# Setup environment ----
library(parallel)
library(data.table)
library(bigKRLS)
library(randomForest)
library(e1071)
library(glmnet)
library(cowplot)
library(ggplot2)
library(ggridges)
library(extrafont)
library(Cairo)
source("R/functions.R")
suppressMessages(loadfonts())
theme_set(theme_minimal() + theme(
  plot.title = element_text(size = 14, hjust = .5),
  text = element_text(family = "CM Roman"),
  panel.border = element_rect(fill = NA, color = "gray70"),
  legend.position = "none"))
results_path <- "results/replication/"
if (!dir.exists(results_path)) {
  dir.create(results_path, recursive = TRUE)
}
if (!dir.exists("figures")) {
  dir.create("figures")
}
if (!dir.exists("tables")) {
  dir.create("tables")
}
load("data/data.RData")
senate_data_lag <- copy(senate_data)
house_data_lag <- copy(house_data)

# Main paper ----
#     Senate data ----
senate_data <- senate_data[jungle == 0 & weird_race == 0 & special == 0 &
    !is.na(ideological_distance) & !is.na(dem_vote) & !is.na(rep_vote) &
    !is.na(real_dem_expenditure_advantage_w_outside)]
sen_tail_lims <- senate_data[,
  quantile(real_dem_expenditure_advantage_w_outside, c(.05, .95),
    na.rm = TRUE)]
sen_X <- data.matrix(senate_data[, .(
  dem_spend_adv = real_dem_expenditure_advantage_w_outside *
    as.numeric(
      real_dem_expenditure_advantage_w_outside >= sen_tail_lims[1] &
      real_dem_expenditure_advantage_w_outside <= sen_tail_lims[2]),
  log_total_spending = log10(real_total_spending_w_outside),
  dem_presvote_advantage,
  adj_dem_presvote_advantage,
  log_dem_presvote = log10(dem_presvote),
  log_rep_presvote = log10(rep_presvote),
  ideological_distance,
  dem_cfscore, rep_cfscore, open_seat, dem_incumbent,
  log_voting_eligible_population = log10(voting_eligible_population),
  y80, y82, y84, y86, y88,
  y90, y92, y94, y96, y98,
  y00, y02, y04, y06, y08,
  y10, y12, y14, y16, y18,
  bottom_tail =
    as.numeric(real_dem_expenditure_advantage_w_outside < sen_tail_lims[1]),
  top_tail =
    as.numeric(real_dem_expenditure_advantage_w_outside > sen_tail_lims[2]),
  bottom_tail_DSA = real_dem_expenditure_advantage_w_outside *
    as.numeric(real_dem_expenditure_advantage_w_outside < sen_tail_lims[1]),
  top_tail_DSA = real_dem_expenditure_advantage_w_outside *
    as.numeric(real_dem_expenditure_advantage_w_outside > sen_tail_lims[2]),
  middle_90 =
    as.numeric(
      real_dem_expenditure_advantage_w_outside <= sen_tail_lims[2] &
      real_dem_expenditure_advantage_w_outside >= sen_tail_lims[1]))])
sen_y <- senate_data[, cbind(dem_vote_share)]

#       data for randomForest ----
sen_X_rf <- data.matrix(senate_data[, .(
  dem_spend_adv = real_dem_expenditure_advantage_w_outside,
  log_total_spending = log10(real_total_spending_w_outside),
  dem_presvote_advantage,
  adj_dem_presvote_advantage,
  log_dem_presvote = log10(dem_presvote),
  log_rep_presvote = log10(rep_presvote),
  ideological_distance,
  dem_cfscore, rep_cfscore, open_seat, dem_incumbent,
  log_voting_eligible_population = log10(voting_eligible_population),
  y80, y82, y84, y86, y88,
  y90, y92, y94, y96, y98,
  y00, y02, y04, y06, y08,
  y10, y12, y14, y16, y18)])
#       dichotomous outcome variable ----
sen_y_01 <- senate_data[, cbind(as.numeric(dem_vote_share > .5))]

#       dataset including lagged outcome ----
senate_data_lag <- merge(senate_data_lag[special == 0],
  senate_data_lag[special == 0, .(year = year + 6, stabb, class,
    lag_outcome = dem_vote / (rep_vote + dem_vote),
    lag_DSA = real_dem_expenditure_advantage_w_outside)],
  by = c("year", "stabb", "class"), all.x = TRUE)
senate_data_lag <- senate_data_lag[
  open_seat == 0 &
    jungle == 0 & weird_race == 0 & special == 0 &
    !is.na(ideological_distance) & !is.na(dem_vote) & !is.na(rep_vote) &
    !is.na(real_dem_expenditure_advantage_w_outside) &
    !is.na(lag_outcome)]
sen_X_lag <- data.matrix(senate_data_lag[, .(
  dem_spend_adv = real_dem_expenditure_advantage_w_outside *
    as.numeric(
      real_dem_expenditure_advantage_w_outside >= sen_tail_lims[1] &
      real_dem_expenditure_advantage_w_outside <= sen_tail_lims[2]),
  log_total_spending = log10(real_total_spending_w_outside),
  dem_presvote_advantage,
  adj_dem_presvote_advantage,
  log_dem_presvote = log10(dem_presvote),
  log_rep_presvote = log10(rep_presvote),
  ideological_distance,
  dem_cfscore, rep_cfscore,
  dem_incumbent,
  log_voting_eligible_population = log10(voting_eligible_population),
  y86, y88,
  y90, y92, y94, y96, y98,
  y00, y02, y04, y06, y08,
  y10, y12, y14, y16, y18,
  bottom_tail =
    as.numeric(real_dem_expenditure_advantage_w_outside < sen_tail_lims[1]),
  top_tail =
    as.numeric(real_dem_expenditure_advantage_w_outside > sen_tail_lims[2]),
  bottom_tail_DSA = real_dem_expenditure_advantage_w_outside *
    as.numeric(real_dem_expenditure_advantage_w_outside < sen_tail_lims[1]),
  top_tail_DSA = real_dem_expenditure_advantage_w_outside *
    as.numeric(real_dem_expenditure_advantage_w_outside > sen_tail_lims[2]),
  middle_90 =
    as.numeric(
      real_dem_expenditure_advantage_w_outside <= sen_tail_lims[2] &
      real_dem_expenditure_advantage_w_outside >= sen_tail_lims[1]),
  lag_outcome = lag_outcome)])
sen_y_lag <- senate_data_lag[, cbind(dem_vote_share)]

#       for numerical marginal effects ----
xS <- sen_X[, 1]
dxS <- setstep(xS)
Xp1S <- make_plus_dx_dataset(sen_X, dxS, "sen")
Xm1S <- make_plus_dx_dataset(sen_X, -dxS, "sen")

#       for CV ----
set.seed(1702567564)
sen_folds <- make_train_and_test_sets(length(sen_y), 5)

#     Senate KRLS ----
if (!file.exists(paste0(results_path, "senate_krls_model"))) {
  senate_krls_model <- bigKRLS(X = sen_X, y = sen_y, noisy = TRUE)
  senate_krls_results <- as.data.table(
    cbind(sen_X, senate_krls_model$derivatives))
  setnames(senate_krls_results,
    c(colnames(sen_X), paste0("d_", colnames(sen_X))))
  senate_krls_cv <- cv_bigKRLS(sen_folds, sen_X, sen_y)
  save.bigKRLS(senate_krls_model,
    model_subfolder_name = paste0(results_path, "senate_krls_model"))
  save(senate_krls_results,
    file = paste0(results_path, "senate_krls_results.RData"))
  save(senate_krls_cv, file = paste0(results_path, "senate_krls_cv.RData"))
} else {
  load.bigKRLS(paste0(results_path, "senate_krls_model"),
    newname = "senate_krls_model")
  load(paste0(results_path, "senate_krls_results.RData"))
  load(paste0(results_path, "senate_krls_cv.RData"))
}

#     Senate KRLS counterfactuals ----
if (!file.exists(paste0(results_path, "senate_krls_counterfactuals.RData"))) {
  set.seed(2009899037)
  qs <- c(list("actual", "zero"), as.list(seq(.95, .05, -.01)))
  seeds <- sample.int(.Machine$integer.max, length(qs))
  senate_krls_counterfactuals <- lapply(seq_along(qs), function(i) {
    type_or_quantile <- qs[[i]]
    set.seed(seeds[[i]])
    cat("\n", type_or_quantile)
    make_counterfactuals(type_or_quantile, senate_krls_model, senate_data)
  })
  save(senate_krls_counterfactuals,
    file = paste0(results_path, "senate_krls_counterfactuals.RData"))
} else {
  load(paste0(results_path, "senate_krls_counterfactuals.RData"))
}

#     Senate figure data ----
if (!file.exists(paste0(results_path, "gg_senate_min_max_data.RData"))) {
  available_seats <- senate_data[ , .(available_seats = .N), year]
  senate_krls_counterfactuals <- rbindlist(senate_krls_counterfactuals)
  senate_krls_counterfactuals <- merge(senate_krls_counterfactuals,
    unincluded_senate_seat_totals[, .(year, D, R)], by = "year")
  gg_senate_min_max_data <- merge(
    senate_krls_counterfactuals[type_or_quantile == "0.95",
      .(boot, year, n_dems_max_dem = n_dems + D)],
    senate_krls_counterfactuals[type_or_quantile == "0.05",
      .(boot, year, n_dems_max_rep = n_dems + D)],
    by = c("boot", "year"))
  gg_senate_min_max_data[, seat_swing := n_dems_max_dem - n_dems_max_rep]
  gg_senate_min_max_data[, n_reps_max_dem := 100 - n_dems_max_dem]
  gg_senate_min_max_data[, n_reps_max_rep := 100 - n_dems_max_rep]
  gg_senate_zero_money_data <- merge(
    senate_krls_counterfactuals[type_or_quantile == "actual",
      .(boot, year, n_dems_actual = n_dems + D)],
    senate_krls_counterfactuals[type_or_quantile == "zero",
      .(boot, year, n_dems_zero = n_dems + D)],
    by = c("boot", "year"))
  gg_senate_zero_money_data[, seat_swing := n_dems_actual - n_dems_zero]
  gg_senate_zero_money_data[, n_reps_actual := 100 - n_dems_actual]
  gg_senate_zero_money_data[, n_reps_zero := 100 - n_dems_zero]
  save(gg_senate_min_max_data,
    file = paste0(results_path, "gg_senate_min_max_data.RData"))
  save(gg_senate_zero_money_data,
    file = paste0(results_path, "gg_senate_zero_money_data.RData"))
} else {
  load(paste0(results_path, "gg_senate_min_max_data.RData"))
  load(paste0(results_path, "gg_senate_zero_money_data.RData"))
}

#     House data ----
house_data <- house_data[!is.na(real_dem_expenditure_advantage) &
    !is.na(ideological_distance) & !is.na(quality_challenger) &
    !is.na(contest_type) & unopposed == 0 & thirdother == 0]
hou_tail_lims <- house_data[,
  quantile(real_dem_expenditure_advantage_w_outside, c(.05, .95), na.rm = TRUE)]
hou_X <- data.matrix(house_data[, .(
  dem_spend_adv = real_dem_expenditure_advantage_w_outside *
    as.numeric(
      real_dem_expenditure_advantage_w_outside >= hou_tail_lims[1] &
        real_dem_expenditure_advantage_w_outside <= hou_tail_lims[2]),
  log_total_spending = log10(real_total_spending_w_outside),
  adj_dem_presvote_advantage, ideological_distance, dem_cfscore, rep_cfscore,
  # house only
  dem_inc_lq_chal = as.numeric(contest_type == 1),
  dem_inc_hq_chal = as.numeric(contest_type == 2),
  rep_inc_lq_chal = as.numeric(contest_type == 3),
  rep_inc_hq_chal = as.numeric(contest_type == 4),
  open_both_hi = as.numeric(contest_type == 5),
  open_hi_dem_lo_rep = as.numeric(contest_type == 6),
  open_lo_dem_hi_rep = as.numeric(contest_type == 7),
  y80, y82, y84, y86, y88,
  y90, y92, y94, y96, y98,
  y00, y02, y04, y06, y08,
  y10, y12, y14, y16, y18,
  bottom_tail =
    as.numeric(real_dem_expenditure_advantage_w_outside < hou_tail_lims[1]),
  top_tail =
    as.numeric(real_dem_expenditure_advantage_w_outside > hou_tail_lims[2]),
  bottom_tail_DSA = real_dem_expenditure_advantage_w_outside *
    as.numeric(real_dem_expenditure_advantage_w_outside < hou_tail_lims[1]),
  top_tail_DSA = real_dem_expenditure_advantage_w_outside *
    as.numeric(real_dem_expenditure_advantage_w_outside > hou_tail_lims[2]),
  middle_90 =
    as.numeric(
      real_dem_expenditure_advantage_w_outside <= hou_tail_lims[2] &
        real_dem_expenditure_advantage_w_outside >= hou_tail_lims[1]))])
hou_y <- house_data[, cbind(dem_vote_share)]
#       data for randomForest ----
hou_X_rf <- data.matrix(house_data[, .(
  dem_spend_adv = real_dem_expenditure_advantage_w_outside,
  log_total_spending = log10(real_total_spending_w_outside),
  adj_dem_presvote_advantage, ideological_distance, dem_cfscore, rep_cfscore,
  dem_inc_lq_chal = as.numeric(contest_type == 1),
  dem_inc_hq_chal = as.numeric(contest_type == 2),
  rep_inc_lq_chal = as.numeric(contest_type == 3),
  rep_inc_hq_chal = as.numeric(contest_type == 4),
  open_both_hi = as.numeric(contest_type == 5),
  open_hi_dem_lo_rep = as.numeric(contest_type == 6),
  open_lo_dem_hi_rep = as.numeric(contest_type == 7),
  y80, y82, y84, y86, y88,
  y90, y92, y94, y96, y98,
  y00, y02, y04, y06, y08,
  y10, y12, y14, y16, y18
)])
#       dichotomous outcome variable ----
hou_y_01 <- house_data[, cbind(as.numeric(dem_vote_share> .5))]
#       data with lagged outcomes ----
house_data_lag <- merge(house_data_lag[redistricted == 0],
  house_data_lag[, .(year = year + 2, district,
    lag_outcome = dem_vote_share,
    lag_DSA = real_dem_expenditure_advantage_w_outside)],
  by = c("year", "district"), all.x = TRUE)
house_data_lag <- house_data_lag[!year %in% c(1980, 1982, 1992, 2002, 2012)]
house_data_lag <- house_data_lag[!is.na(real_dem_expenditure_advantage) &
    !is.na(ideological_distance) & !is.na(quality_challenger) &
    !is.na(contest_type) & unopposed == 0 & thirdother == 0 &
    !is.na(lag_outcome)]
hou_X_lag <- data.matrix(house_data_lag[, .(
  dem_spend_adv = real_dem_expenditure_advantage_w_outside *
    as.numeric(
      real_dem_expenditure_advantage_w_outside >= hou_tail_lims[1] &
      real_dem_expenditure_advantage_w_outside <= hou_tail_lims[2]),
  log_total_spending = log10(real_total_spending_w_outside),
  adj_dem_presvote_advantage, ideological_distance, dem_cfscore, rep_cfscore,
  # house only
  dem_inc_lq_chal = as.numeric(contest_type == 1),
  dem_inc_hq_chal = as.numeric(contest_type == 2),
  rep_inc_lq_chal = as.numeric(contest_type == 3),
  rep_inc_hq_chal = as.numeric(contest_type == 4),
  y84, y86,
  y88, y90,
  y94,
  y96, y98,
  y00,
  y04, y06,
  y08, y10,
  y14,
  y16, y18,
  bottom_tail =
    as.numeric(real_dem_expenditure_advantage_w_outside < hou_tail_lims[1]),
  top_tail =
    as.numeric(real_dem_expenditure_advantage_w_outside > hou_tail_lims[2]),
  bottom_tail_DSA = real_dem_expenditure_advantage_w_outside *
    as.numeric(real_dem_expenditure_advantage_w_outside < hou_tail_lims[1]),
  top_tail_DSA = real_dem_expenditure_advantage_w_outside *
    as.numeric(real_dem_expenditure_advantage_w_outside > hou_tail_lims[2]),
  middle_90 =
    as.numeric(
      real_dem_expenditure_advantage_w_outside <= hou_tail_lims[2] &
        real_dem_expenditure_advantage_w_outside >= hou_tail_lims[1]),
  lag_outcome = lag_outcome)])
hou_y_lag <- house_data_lag[, cbind(dem_vote_share)]

#       setup data for marginal effects ----
xH <- hou_X[, 1]
dxH <- setstep(xH)
Xp1H <- make_plus_dx_dataset(hou_X, dxH, "hou")
Xm1H <- make_plus_dx_dataset(hou_X, -dxH, "hou")
#       for CV ----
set.seed(1546625538)
hou_folds <- make_train_and_test_sets(length(hou_y), 5)

#     House KRLS ----
if (!file.exists(paste0(results_path, "house_krls_results.RData"))) {
  house_krls_model <- bigKRLS::bigKRLS(X = hou_X, y = hou_y, noisy = TRUE)
  house_krls_results <- as.data.table(cbind(hou_X, house_krls_model$derivatives))
  setnames(house_krls_results, c(colnames(hou_X), paste0("d_", colnames(hou_X))))
  house_krls_cv <- cv_bigKRLS(hou_folds, hou_X, hou_y)
  save.bigKRLS(house_krls_model,
    model_subfolder_name = paste0(results_path, "house_krls_model"))
  save(house_krls_results,
    file = paste0(results_path, "house_krls_results.RData"))
  save(house_krls_cv, file = paste0(results_path, "house_krls_cv.RData"))
} else {
  load.bigKRLS(paste0(results_path, "house_krls_model"),
    newname = "house_krls_model")
  load(paste0(results_path, "house_krls_results.RData"))
  load(paste0(results_path, "house_krls_cv.RData"))
}

#     House KRLS counterfactuals ----
if (!file.exists(paste0(results_path, "house_krls_counterfactuals.RData"))) {
  RNGkind("L'Ecuyer")
  qs <- c(list("actual", "zero"), as.list(seq(.95, .05, -.01)))
  set.seed(1961685937)
  seeds <- sample.int(.Machine$integer.max, length(qs))
  house_krls_counterfactuals <- mclapply(seq_along(qs), function(i) {
    type_or_quantile <- qs[[i]]
    set.seed(seeds[[i]])
    cat("\n", format(Sys.time()), type_or_quantile)
    make_counterfactuals(type_or_quantile, house_krls_model, house_data)
  }, mc.cores = 6)
  save(house_krls_counterfactuals,
    file = paste0(results_path, "house_krls_counterfactuals.RData"))
} else {
  load(paste0(results_path, "house_krls_counterfactuals.RData"))
}

#     House figure data ----
if (!file.exists(paste0(results_path, "gg_house_min_max_data.RData"))) {
  available_seats <- house_data[ , .(available_seats = .N), year]
  house_counterfactuals <- rbindlist(house_krls_counterfactuals)
  house_counterfactuals <- merge(house_counterfactuals,
    unincluded_house_seat_totals[, .(year, D, R)], by = "year")
  gg_house_min_max_data <- merge(
    house_counterfactuals[type_or_quantile == "0.95",
      .(boot, year, n_dems_max_dem = n_dems + D)],
    house_counterfactuals[type_or_quantile == "0.05",
      .(boot, year, n_dems_max_rep = n_dems + D)],
    by = c("boot", "year"))
  gg_house_min_max_data[, seat_swing := n_dems_max_dem - n_dems_max_rep]
  gg_house_min_max_data[, n_reps_max_dem := 435 - n_dems_max_dem]
  gg_house_min_max_data[, n_reps_max_rep := 435 - n_dems_max_rep]
  gg_house_zero_money_data <- merge(
    house_counterfactuals[type_or_quantile == "actual",
      .(boot, year, n_dems_actual = n_dems + D)],
    house_counterfactuals[type_or_quantile == "zero",
      .(boot, year, n_dems_zero = n_dems + D)],
    by = c("boot", "year"))
  gg_house_zero_money_data[, seat_swing := n_dems_actual - n_dems_zero]
  gg_house_zero_money_data[, n_reps_actual := 435 - n_dems_actual]
  gg_house_zero_money_data[, n_reps_zero := 435 - n_dems_zero]
  save(gg_house_min_max_data,
    file = paste0(results_path, "gg_house_min_max_data.RData"))
  save(gg_house_zero_money_data,
    file = paste0(results_path, "gg_house_zero_money_data.RData"))
} else {
  load(file = paste0(results_path, "gg_house_min_max_data.RData"))
  load(file = paste0(results_path, "gg_house_zero_money_data.RData"))
}

#     Make Figure 1 ----
house_krls_results[, `:=`(
  group = factor(dplyr::case_when(
    dem_inc_lq_chal + dem_inc_hq_chal == 1 ~ "Dem. Incumbent",
    rep_inc_lq_chal + rep_inc_hq_chal == 1 ~ "Rep. Incumbent",
    TRUE ~ "All"
  ), levels = c("All",
    "Rep. Incumbent", "Dem. Incumbent"))
)]
senate_krls_results[, `:=`(
  group = factor(dplyr::case_when(
    dem_incumbent == 1 ~ "Dem. Incumbent",
    dem_incumbent + open_seat == 0 ~ "Rep. Incumbent",
    TRUE ~ "All"
  ), levels = c("All",
    "Rep. Incumbent", "Dem. Incumbent"))
)]

fig_house_nonmonotonic <-
  ggplot(house_krls_results[bottom_tail == 0 & top_tail == 0],
    aes(dem_spend_adv, d_dem_spend_adv)) +
  geom_hline(yintercept = 0, lty = 3, color = "gray") +
  geom_point(alpha = .05, shape = 21, fill = "black", color = "white",
    stroke = 0) +
  geom_smooth(data = house_krls_results[bottom_tail == 0 & top_tail == 0 &
      group != "All"],
    aes(group = group, linetype = group),
    method = "loess", se = FALSE) +
  geom_smooth(method = "loess") +
  scale_linetype_manual(values = c(
    "All" = "solid",
    "Rep. Incumbent" = "dotdash",
    "Dem. Incumbent" = "dashed"), drop = FALSE) +
  scale_x_continuous(breaks = seq(-1.5, 1.5, .5)) +
  scale_y_continuous(breaks = c(0, .03, .06),
    labels = c("0%", "3%", "6%")) +
  xlab(" ") + ylab("Marginal Effect of Spending\non Democratic Vote Share") +
  ggtitle("House") +
  coord_cartesian(expand = FALSE, xlim = c(-1.75, 1.75), ylim = c(-.015, .06)) +
  theme(legend.position = "bottom", legend.title = element_blank())
fig_senate_nonmonotonic <-
  ggplot(senate_krls_results[bottom_tail == 0 & top_tail == 0],
    aes(dem_spend_adv, d_dem_spend_adv)) +
  geom_hline(yintercept = 0, lty = 3, color = "gray") +
  geom_point(alpha = .2, shape = 21, fill = "black", color = "white",
    stroke = 0) +
  geom_smooth(data = senate_krls_results[bottom_tail == 0 & top_tail == 0 &
      group != "All"], aes(group = group, linetype = group),
    method = "loess", se = FALSE) +
  geom_smooth(method = "loess") +
  scale_linetype_manual(values = c(
    "All" = "solid",
    "Rep. Incumbent" = "dotdash",
    "Dem. Incumbent" = "dashed")) +
  scale_x_continuous(breaks = seq(-10, 10, 5)) +
  scale_y_continuous(breaks = c(0, .005, .01, .015),
    labels = c("0%", "0.5%", "1%", "1.5%")) +
  xlab(" ") + ylab(" ") + ggtitle("Senate") +
  coord_cartesian(expand = FALSE, xlim = c(-11, 11), ylim = c(-.015/4, .015))
h <- .125
h2 <- .03
fig_nonmonotonic <- cowplot::ggdraw() +
  cowplot::draw_plot(fig_house_nonmonotonic + theme(
    legend.position = "none",
    plot.title = element_text(size = 12),
    text = element_text(family = "CM Roman")),
    x = 0, y = h,  height = 1 - 2 * h, width = .5) +
  cowplot::draw_plot(fig_senate_nonmonotonic + theme(
    plot.title = element_text(size = 12),
    text = element_text(family = "CM Roman")),
    x = .5, y = h, height = 1 - 2 * h, width = .5 - h2) +
  cowplot::draw_label(
    paste0("Money Is Most Effective with Spending Near Parity"),
    size = 18, x = .5, y = 1 - h / 2,
    fontface = "bold",
    fontfamily = "CM Roman") +
  cowplot::draw_label(paste0(
    "Democratic Candidate's Spending Advantage\n(in millions of ",
    "inflation-adjusted US$)"),
    size = 12, x = .5, y = .75 * h+.025,
    fontfamily = "CM Roman") +
  cowplot::draw_plot(
    get_legend(fig_house_nonmonotonic + theme(
      plot.title = element_text(size = 12),
      text = element_text(family = "CM Roman"))),
    x = 0, y = 0, height = .05, width = 1)

pdf(paste0("figures/fig-nonmonotonic-", Sys.Date(), ".pdf"),
  width = 6.5, height = 4)
print(fig_nonmonotonic)
dev.off()
embed_fonts(paste0("figures/fig-nonmonotonic-", Sys.Date(), ".pdf"))


#     Make Figure 2 ----
fig_house_min_max <- ggplot(
  gg_house_min_max_data,
  aes(y = year, group = year, fill = ..x..)) +
  geom_vline(xintercept = 218, lty = 3, color = "gray") +
  geom_density_ridges(
    aes(x = n_dems_max_dem),
    bandwidth = 4.35,
    scale = 1, size = 0.1, rel_min_height = 0.01,
    color = "gray92",
    #alpha = .25,
    fill = "#0276FD") +
  geom_density_ridges(
    aes(x = n_dems_max_rep),
    bandwidth = 4.35,
    scale = 1, size = 0.1, rel_min_height = 0.01,
    color = "gray92",
    alpha = .5,
    fill = "#ef2f4f") +
  scale_x_continuous(limits = c(50, 385), expand = c(0, 0),
    breaks = c(70, 145, 218, 290, 360)) +
  scale_y_reverse(breaks=seq(2016, 1980, -4), expand = c(0, 0.2)) +
  xlab("Democratic House Seats") +
  theme_minimal() +
  theme(panel.border = element_rect(color = "gray92", fill = NA),
    legend.position = "none",
    plot.title = element_text(hjust = .5)) +
  ggtitle("") +
  ylab("")

fig_senate_min_max <- ggplot(
  gg_senate_min_max_data,
  aes(y = year, group = year, fill = ..x..)) +
  geom_vline(xintercept = 50, lty = 3, color = "gray") +
  geom_density_ridges(
    aes(x = n_dems_max_dem),
    bandwidth = 1,
    scale = 1, size = 0.1, rel_min_height = 0.01,
    color = "gray92",
    #alpha = .25,
    fill = "#0276FD") +
  geom_density_ridges(
    aes(x = n_dems_max_rep),
    bandwidth = 1,
    scale = 1, size = 0.1, rel_min_height = 0.01,
    color = "gray92",
    alpha = .5,
    fill = "#ef2f4f") +
  scale_x_continuous(limits = c(25, 75), expand = c(0, 0),
    breaks = seq(30, 70, 10)) +
  scale_y_reverse(breaks=seq(2016, 1980, -4), expand = c(0, 0.2)) +
  xlab("Democratic Senate Seats") +
  theme_minimal() +
  theme(panel.border = element_rect(color = "gray92", fill = NA),
    legend.position = "none",
    plot.title = element_text(hjust = .5)) +
  ggtitle("") +
  ylab("")

h <- .15
h2 <- .03
fig_min_max <-
  cowplot::ggdraw() +
    cowplot::draw_plot(fig_house_min_max + ggtitle("House") +
        xlab("") + theme(
          plot.title = element_text(size = 12),
          text = element_text(family = "CM Roman")) +
        scale_y_reverse(breaks = seq(2016, 1980, -6),
          minor_breaks = seq(2016, 1980, -2), expand = c(0, .3)) +
        theme(panel.grid.minor = element_line(size = .1),
          panel.grid.major.y = element_line(size = .1)),
      x = 0, y = .35 * h,  height = 1 - 1 * h, width = .5) +
    cowplot::draw_plot(fig_senate_min_max + ggtitle("Senate") +
        xlab("") + theme(
          plot.title = element_text(size = 12),
          text = element_text(family = "CM Roman")) +
        scale_y_reverse(breaks = seq(2016, 1980, -6),
          minor_breaks = seq(2016, 1980, -2), expand = c(0, .3)) +
        theme(panel.grid.minor = element_line(size = .1),
          panel.grid.major.y = element_line(size = .1)),
      x = .5, y = .35 * h, height = 1 - 1 * h, width = .5 - h2) +
    cowplot::draw_label(
      paste0("Congressional Control Depends on Spending"),
      size = 18, x = .5, y = 1 - .4 * h, fontface = "bold",
      fontfamily = "CM Roman") +
    cowplot::draw_label(paste0(
      "Number of Seats Held by Democrats Under Hypothetical Spending Profiles"),
      size = 12, x = .5, y = .4 * h,
      fontfamily = "CM Roman")

pdf(paste0("figures/fig-min-max-", Sys.Date(), ".pdf"),
  width = 6.5, height = 4)
print(fig_min_max)
dev.off()
embed_fonts(paste0("figures/fig-min-max-", Sys.Date(), ".pdf"))

#     Make Figure 3 ----
house_sig_years <- gg_house_zero_money_data[,
  .(year, V1 = jitter(n_dems_zero - n_dems_actual, factor = .001))
  ][, .(mean(V1), quantile(V1, .025), quantile(V1, .975)), year][
    which(sign(V2) == sign(V3) | V2 == 0 | V3 == 0), year]
senate_sig_years <- gg_senate_zero_money_data[,
  .(year, V1 = jitter(n_dems_zero - n_dems_actual, factor = .001))
  ][, .(mean(V1), quantile(V1, .025), quantile(V1, .975)), year][
    which(sign(V2) == sign(V3) | V2 == 0 | V3 == 0), year]
gg_house_zero_money_data[,
  sig := as.numeric(year %in% house_sig_years)]
gg_senate_zero_money_data[,
  sig := as.numeric(year %in% senate_sig_years)]
fig_house_zero_money <- ggplot(
  gg_house_zero_money_data,
  aes(y = year, group = year, fill = ..x..)) +
  geom_vline(xintercept = 0, lty = 3, color = "gray") +
  geom_density_ridges(
    aes(x = n_dems_zero - n_dems_actual, fill = factor(sig)),
    bandwidth = 2.35,
    scale = 1, size = 0.1, rel_min_height = 0.05,
    color = "gray92") +
  scale_y_reverse(breaks=seq(2016, 1980, -4), expand = c(0, 0.2)) +
  xlab("Change in Democratic House Seats") +
  theme_minimal() +
  theme(panel.border = element_rect(color = "gray92", fill = NA),
    legend.position = "none",
    plot.title = element_text(hjust = .5)) +
  ggtitle("") +
  scale_fill_manual(values = c("0" = "gray72", "1" = "gray42")) +
  ylab("")

fig_senate_zero_money <- ggplot(
  gg_senate_zero_money_data,
  aes(y = year, group = year, fill = ..x..)) +
  geom_vline(xintercept = 0, lty = 3, color = "gray") +
  geom_density_ridges(
    aes(x = n_dems_zero - n_dems_actual, fill = factor(sig)),
    bandwidth = 1,
    scale = 1, size = 0.1, rel_min_height = 0.05,
    color = "gray92") +
  scale_y_reverse(breaks=seq(2016, 1980, -4), expand = c(0, 0.2)) +
  xlab("Change in Democratic Senate Seats") +
  theme_minimal() +
  theme(panel.border = element_rect(color = "gray92", fill = NA),
    legend.position = "none",
    plot.title = element_text(hjust = .5)) +
  ggtitle("") +
  scale_fill_manual(values = c("0" = "gray72", "1" = "gray42")) +
  ylab("")

h <- .15
h2 <- .03
fig_zero_money <- cowplot::ggdraw() +
    cowplot::draw_plot(fig_house_zero_money + ggtitle("House") +
        xlab("") + theme(
          plot.title = element_text(size = 12),
          text = element_text(family = "CM Roman")) +
        scale_y_reverse(expand = c(0, .3)),
      x = 0, y = .35 * h,  height = 1 - 1.5 * h, width = .5) +
    cowplot::draw_plot(fig_senate_zero_money + ggtitle("Senate") +
        xlab("") + theme(
          plot.title = element_text(size = 12),
          text = element_text(family = "CM Roman")) +
        scale_y_reverse(expand = c(0, .3)),
      x = .5, y = .35 * h, height = 1 - 1.5 * h, width = .5 - h2) +
    cowplot::draw_label(
      paste0("The Effect of Removing Money on Democratic Seats"),
      size = 18, x = .5, y = 1 - .8 * h, fontface = "bold",
      fontfamily = "CM Roman") +
    cowplot::draw_label(paste0(
      "Additional Seats Held by Democrats under Zero Spending"),
      size = 12, x = .5, y = .25 * h,
      fontfamily = "CM Roman")

pdf(paste0("figures/fig-zero-money-", Sys.Date(), ".pdf"),
  width = 6.5, height = 4)
print(fig_zero_money)
dev.off()
embed_fonts(paste0("figures/fig-zero-money-", Sys.Date(), ".pdf"))

#     Make Figure 4 ----
house_gap_data <- house_data[
  real_dem_expenditure_advantage_w_outside >= hou_tail_lims[1] &
    real_dem_expenditure_advantage_w_outside <= hou_tail_lims[2],
  .(gap = mean(real_dem_expenditure_advantage_w_outside[
    real_dem_expenditure_advantage_w_outside > 0]) -
      mean(real_dem_expenditure_advantage_w_outside[
        real_dem_expenditure_advantage_w_outside < 0])),
  year]
senate_gap_data <- senate_data[
  real_dem_expenditure_advantage_w_outside >= sen_tail_lims[1] &
    real_dem_expenditure_advantage_w_outside <= sen_tail_lims[2],
  .(gap = mean(real_dem_expenditure_advantage_w_outside[
    real_dem_expenditure_advantage_w_outside > 0]) -
      mean(real_dem_expenditure_advantage_w_outside[
        real_dem_expenditure_advantage_w_outside < 0])),
  year]

fig_house_dmr <-
  ggplot(house_krls_results[bottom_tail == 0 & top_tail == 0],
    aes(10^log_total_spending, d_dem_spend_adv)) +
  geom_hline(yintercept = 0, lty = 3, color = "gray") +
  geom_point(alpha = .05, shape = 21, fill = "black", color = "white",
    stroke = 0) +
  geom_smooth(method = "loess") +
  scale_x_log10(breaks = c(1e5, 1e6, 1e7), labels = c("$100k", "$1M", "$10M")) +
  scale_y_continuous(breaks = c(0, .03, .06),
    labels = c("0%", "3%", "6%")) +
  xlab(" ") + ylab("Marginal Effect of Spending\non Democratic Vote Share") +
  ggtitle(" ") +
  coord_cartesian(expand = FALSE, xlim = c(1e5-1e3, 1e7+1e5), ylim = c(-.015, .07)) +
  theme(legend.position = "bottom", legend.title = element_blank())

fig_senate_dmr <-
  ggplot(senate_krls_results[bottom_tail == 0 & top_tail == 0],
    aes(10^log_total_spending, d_dem_spend_adv)) +
  geom_hline(yintercept = 0, lty = 3, color = "gray") +
  geom_point(alpha = .2, shape = 21, fill = "black", color = "white",
    stroke = 0) +
  geom_smooth(method = "loess") +
  scale_x_log10(breaks = c(1e6, 1e7, 1e8), labels = c("$1M", "$10M", "$100M")) +
  scale_y_continuous(breaks = c(0, .005, .01, .015),
    labels = c("0%", "0.5%", "1%", "1.5%")) +
  xlab(" ") + ylab(" ") + ggtitle(" ") +
  coord_cartesian(expand = FALSE, xlim = c(1e6-1e4, 1e8+1e6), ylim = c(-.015/(4.667), .015)) +
  theme(legend.position = "bottom", legend.title = element_blank())

house_gap_data <- house_data[
  real_dem_expenditure_advantage_w_outside >= hou_tail_lims[1] &
    real_dem_expenditure_advantage_w_outside <= hou_tail_lims[2],
  .(total = mean(real_total_spending_w_outside) / 1e6,
    gap = mean(abs(real_dem_expenditure_advantage_w_outside))
  ), year]

senate_gap_data <- senate_data[
  real_dem_expenditure_advantage_w_outside >= sen_tail_lims[1] &
    real_dem_expenditure_advantage_w_outside <= sen_tail_lims[2],
  .(total = mean(real_total_spending_w_outside) / 1e6,
    gap = mean(abs(real_dem_expenditure_advantage_w_outside))
  ), year]

fig_house_gap <- ggplot(house_gap_data, aes(year, gap)) +
  geom_point() +
  geom_smooth(method = "lm") +
  ylab("Mean Expenditure Advantage\n(Real Millions US$)") +
  xlab("") +
  scale_y_continuous(limits = c(0, 1.1),
    breaks = c(0, .5, 1), labels = c(" $0M", "$0.5M", " $1M")) +
  coord_cartesian(xlim = c(1979, 2019), ylim = c(0, 1.1), expand = FALSE) +
  ggtitle("House")

fig_senate_gap <- ggplot(senate_gap_data, aes(year, gap)) +
  geom_point() +
  geom_smooth(method = "lm") +
  ylab(" ") +
  xlab("") +
  scale_y_continuous(limits = c(0, 5.5),
    breaks = c(0, 2, 4), labels = c(" $0M", " $2M", " $4M")) +
  coord_cartesian(xlim = c(1979, 2019), expand = FALSE) +
  ggtitle("Senate")

fig_house_total <- ggplot(house_gap_data, aes(year, total)) +
  geom_point() +
  geom_smooth(method = "lm") +
  ylab("Mean Total Spending\n(Real Millions US$)") +
  # ylab("Millions of Inflation-adjusted US$") +
  xlab(" ") +
  scale_y_continuous(limits = c(0, 3.1),
    breaks = c(0, 1, 2, 3), labels = c("$0M", " $1M", " $2M", " $3M")) +
  coord_cartesian(xlim = c(1979, 2019), expand = FALSE) +
  ggtitle(" ")

fig_senate_total <- ggplot(senate_gap_data, aes(year, total)) +
  geom_point() +
  geom_smooth(method = "lm") +
  ylab(" ") +
  xlab("") +
  scale_y_continuous(limits = c(0, 26),
    breaks = c(0, 10, 20), labels = c("$0M", "$10M", "$20M")) +
  coord_cartesian(xlim = c(1979, 2019), expand = FALSE) +
  ggtitle(" ")

h <- .8 / 3
w <- .03
fig_gap_dmr <- cowplot::ggdraw() +
    cowplot::draw_plot(fig_house_gap,
      x = 0, y = .05 + 2 * h,  height = h, width = .5) +
    cowplot::draw_plot(fig_senate_gap,
      x = .5, y = .05 + 2 * h, height = h, width = .5 - w) +
    cowplot::draw_plot(fig_house_total,
      x = 0, y = .05 + 1 * h,  height = h, width = .5) +
    cowplot::draw_plot(fig_senate_total,
      x = .5, y = .05 + 1 * h, height = h, width = .5 - w) +
    cowplot::draw_plot(fig_house_dmr,
      x = 0, y = .05,  height = h, width = .5) +
    cowplot::draw_plot(fig_senate_dmr,
      x = .5, y = .05, height = h, width = .5 - w) +
    cowplot::draw_label(
      paste0("Vanishing Marginals and Decreasing Returns\nAttenuate the Effectiveness of Spending"),
      size = 18, x = .5, y = .925, fontface = "bold",
      fontfamily = "CM Roman") +
    cowplot::draw_label(paste0(
      "Log of Total Spending"),
      size = 10, x = .5, y = .06,
      fontfamily = "CM Roman")

pdf(paste0("figures/fig-gap-dmr-", Sys.Date(), ".pdf"),
  width = 6.5, height = 9)
print(fig_gap_dmr)
dev.off()
embed_fonts(paste0("figures/fig-gap-dmr-", Sys.Date(), ".pdf"))

# Reported quantities ----
#     p. 6, integrating areas under marginal effects curves ----
loess_house <- house_krls_results[bottom_tail + top_tail == 0,
  loess(d_dem_spend_adv ~ dem_spend_adv)]
x_hou <- integrate(function(x) predict(loess_house, x),
  lower = house_krls_results[bottom_tail + top_tail == 0, min(dem_spend_adv)],
  upper = house_krls_results[bottom_tail + top_tail == 0, max(dem_spend_adv)])
x_hou

loess_senate <- senate_krls_results[bottom_tail + top_tail == 0,
  loess(d_dem_spend_adv ~ dem_spend_adv)]
x_sen <- integrate(function(x) predict(loess_senate, x),
  lower = senate_krls_results[bottom_tail + top_tail == 0, min(dem_spend_adv)],
  upper = senate_krls_results[bottom_tail + top_tail == 0, max(dem_spend_adv)])
x_sen

# marginal seat proportions:
house_data[year == 2018, mean(
  dem_vote / (dem_vote + rep_vote) > .5 - x_hou$value / 2 &
    dem_vote / (dem_vote + rep_vote) < .5 + x_hou$value / 2)]
senate_data[year == 2018, mean(
  dem_vote / (dem_vote + rep_vote) > .5 - x_sen$value / 2 &
    dem_vote / (dem_vote + rep_vote) < .5 + x_sen$value / 2)]

#     p. 7, back-of-the-envelope for linear extrapolation ----
# $100k in 1998 dollars converted to 2016 dollars:
unit_cost <- 100000 / house_data[year == 1998, year_scale[1]]

# from -$1.7M to $1.7M = 3.4e6
overall_added_votes <- 3.4e6 / unit_cost

# upper estimate of additional seats from linear extrapolation
2.17 * overall_added_votes


#     p. 7, scaled-up system-level costs ----
# on average for dems:
house_data[
  real_dem_expenditure_advantage_w_outside > hou_tail_lims[1] &
    real_dem_expenditure_advantage_w_outside < hou_tail_lims[2],
  sum(hou_tail_lims[2] - real_dem_expenditure_advantage_w_outside),
  year][, mean(V1)]
# on average for reps:
house_data[
  real_dem_expenditure_advantage_w_outside > hou_tail_lims[1] &
    real_dem_expenditure_advantage_w_outside < hou_tail_lims[2],
  sum(hou_tail_lims[1] - real_dem_expenditure_advantage_w_outside),
  year][, mean(V1)]
# and in the senate:
# for dems:
senate_data[
  real_dem_expenditure_advantage_w_outside > sen_tail_lims[1] &
    real_dem_expenditure_advantage_w_outside < sen_tail_lims[2],
  sum(sen_tail_lims[2] - real_dem_expenditure_advantage_w_outside),
  year][, mean(V1)]
# for reps:
senate_data[
  real_dem_expenditure_advantage_w_outside > sen_tail_lims[1] &
    real_dem_expenditure_advantage_w_outside < sen_tail_lims[2],
  sum(sen_tail_lims[1] - real_dem_expenditure_advantage_w_outside),
  year][, mean(V1)]
# current total House spending by Dems
house_data[year == 2016 &
    real_dem_expenditure_advantage_w_outside >= sen_tail_lims[1] &
    real_dem_expenditure_advantage_w_outside <= sen_tail_lims[2],
  sum(real_dem_expenditure_w_outside / year_scale / 1e6)]
# current total House spending by Reps
house_data[year == 2016 &
    real_dem_expenditure_advantage_w_outside >= sen_tail_lims[1] &
    real_dem_expenditure_advantage_w_outside <= sen_tail_lims[2],
  sum(real_rep_expenditure_w_outside / year_scale / 1e6),
  year][, mean(V1)]
# current total Senate spending by Dems
senate_data[year == 2016 &
    real_dem_expenditure_advantage_w_outside >= sen_tail_lims[1] &
    real_dem_expenditure_advantage_w_outside <= sen_tail_lims[2],
  sum(real_dem_expenditure_w_outside / year_scale / 1e6)]
# current total Senate spending by Reps
senate_data[year == 2016 &
    real_dem_expenditure_advantage_w_outside >= sen_tail_lims[1] &
    real_dem_expenditure_advantage_w_outside <= sen_tail_lims[2],
  sum(real_rep_expenditure_w_outside / year_scale / 1e6),
  year][, mean(V1)]

#     p. 9, minimum observed spending levels ----
# min total spending, excl. tail cases
round(10 ^ min(hou_X[hou_X[, 34] + hou_X[, 35] == 0, "log_total_spending"]), -3)
round(10 ^ min(sen_X[sen_X[, 33] + sen_X[, 34] == 0, "log_total_spending"]), -3)


# Appendix A ----
round(sen_tail_lims, 1)
round(hou_tail_lims, 1)

# Appendix C ----
#     Ranges and Kurtosis ----
round(quantile(senate_data$real_dem_expenditure_advantage_w_outside,
  c(.05, .95), na.rm = TRUE), 1)
round(range(senate_data$real_dem_expenditure_advantage_w_outside, na.rm = TRUE))
round(e1071::kurtosis(senate_data$real_dem_expenditure_advantage_w_outside,
  na.rm = TRUE, type = 1), 1)

round(quantile(house_data$real_dem_expenditure_advantage_w_outside, c(.05, .95),
  na.rm = TRUE), 1)
round(range(house_data$real_dem_expenditure_advantage_w_outside, na.rm = TRUE))
round(e1071::kurtosis(house_data$real_dem_expenditure_advantage_w_outside,
  na.rm = TRUE, type = 1), 1)

#     House KRLS results for appendix ----
house_krls_sum <- summary(house_krls_model)
house_krls_sum$ttests <- as.data.table(house_krls_sum$ttests)
house_krls_sum$percentiles <- as.data.table(house_krls_sum$percentiles)
hou_index <- rbind(
  c(1, "Democratic Spending Advantage"),
  c(2, "$\\text{log}$(Total Spending)"),
  c(3, "Democratic Presidential Vote Advantage"),
  c(4, "Relative Republican Extremism"),
  c(5, "Democrat Candidate's CF Score"),
  c(6, "Republican Candidate's CF Score"),
  c(7, "Dem. Inc./Low Qual. Chall."),
  c(8, "Dem. Inc./High Qual. Chall."),
  c(9, "Rep Inc./Low Qual. Chall."),
  c(10, "Rep Inc./High Qual. Chall."),
  c(11, "Open Seat/Both High Qual."),
  c(12, "Open Seat/High Qual. Dem/Low Qual. Rep"),
  c(13, "Open Seat/Low Qual. Dem/High Qual. Rep"),
  c(14, "Year = 1980"),
  c(15, "Year = 1982"),
  c(16, "Year = 1984"),
  c(17, "Year = 1986"),
  c(18, "Year = 1988"),
  c(19, "Year = 1990"),
  c(20, "Year = 1992"),
  c(21, "Year = 1994"),
  c(22, "Year = 1996"),
  c(23, "Year = 1998"),
  c(24, "Year = 2000"),
  c(25, "Year = 2002"),
  c(26, "Year = 2004"),
  c(27, "Year = 2006"),
  c(28, "Year = 2008"),
  c(29, "Year = 2010"),
  c(30, "Year = 2012"),
  c(31, "Year = 2014"),
  c(32, "Year = 2016"),
  c(33, "Year = 2018"),
  c(34, "Bottom 5\\%"),
  c(38, "Middle 90\\%"),
  c(35, "Top 5\\%"),
  c(36, "Bottom 5\\% $\\times$ Dem. Spending Adv."),
  c(37, "Top 5\\% $\\times$ Dem. Spending Adv.")
)
sink("tables/table-s1.txt")
for (i in 1:nrow(hou_index)) {
  make_outcome_row(house_krls_sum$ttests[as.numeric(hou_index[i, 1]), ],
    hou_index[i, 2])
}
sink()

#     Senate KRLS results for appendix ----
senate_krls_sum <- summary(senate_krls_model)
senate_krls_sum$ttests <- as.data.table(senate_krls_sum$ttests)
senate_krls_sum$percentiles <- as.data.table(senate_krls_sum$percentiles)
sen_index <- rbind(
  c(1, "Democratic Spending Advantage"),
  c(2, "$\\text{log}$(Total Spending)"),
  c(3, "Democratic Presidential Vote Advantage"),
  c(4, "Adj. Dem. Pres. Vote Advantage"),
  c(5, "$\\text{log}$($n$ Votes for Democratic Presidential Candidate)"),
  c(6, "$\\text{log}$($n$ Votes for Republican Presidential Candidate)"),
  c(7, "Relative Republican Extremism"),
  c(8, "Democrat Candidate's CF Score"),
  c(9, "Republican Candidate's CF Score"),
  c(10, "Open Seat"), c(11, "Democratic Incumbent"),
  c(12, "$\\text{log}$(Voting Eligible Population)"),
  c(13, "Year = 1980"), c(14, "Year = 1982"), c(15, "Year = 1984"),
  c(16, "Year = 1986"), c(17, "Year = 1988"), c(18, "Year = 1990"),
  c(19, "Year = 1992"), c(20, "Year = 1994"), c(21, "Year = 1996"),
  c(22, "Year = 1998"), c(23, "Year = 2000"), c(24, "Year = 2002"),
  c(25, "Year = 2004"), c(26, "Year = 2006"), c(27, "Year = 2008"),
  c(28, "Year = 2010"), c(29, "Year = 2012"), c(30, "Year = 2014"),
  c(31, "Year = 2016"), c(32, "Year = 2018"),
  c(33, "Bottom 5\\%"),
  c(37, "Middle 90\\%"),
  c(34, "Top 5\\%"),
  c(35, "Bottom 5\\% $\\times$ Dem. Spending Adv."),
  c(36, "Top 5\\% $\\times$ Dem. Spending Adv.")
)
sink("tables/table-s2.txt")
for (i in 1:nrow(sen_index)) {
  make_outcome_row(senate_krls_sum$ttests[as.numeric(sen_index[i, 1]), ],
    sen_index[i, 2])
}
sink()

# Appendix D ----
#     Descriptive stats tables ----
house_data_trans <- copy(house_data)
senate_data_trans <- copy(senate_data)
load("data.RData")
sink("tables/table-s3.txt")
house_data[, rbind(
  make_description_row(dem_vote_share, "Democratic Vote Share"),
  make_description_row(real_dem_expenditure_advantage_w_outside,
    "Democratic Expenditure Advantage (Millions of US\\$)"),
  make_description_row(log10(real_total_spending_w_outside),
    "(log) Total Expenditure (Millions of US\\$)"),
  make_description_row(dem_cfscore, "Democrat CF Score"),
  make_description_row(rep_cfscore, "Republican CF Score"),
  make_description_row(ideological_distance, "Relative Republican Extremism"),
  make_description_row(adj_dem_presvote_advantage,
    "Adjusted Democratic Presidential Vote Advantage"),
  make_description_row(open_seat, "Open Seat"),
  make_description_row(dem_incumbent, "Democratic Incumbent"),
  make_description_row(quality_challenger, "Quality Challenger"),
  make_description_row(unopposed, "Unopposed"))]
house_data[, .N]
nrow(hou_X)
sink()
sink("tables/table-s4.txt")
senate_data[special == 0, rbind(
  make_description_row(dem_vote_share, "Democratic Vote Share"),
  make_description_row(real_dem_expenditure_advantage_w_outside,
    "Democratic Expenditure Advantage (Millions of US\\$)"),
  make_description_row(log10(real_total_spending_w_outside),
    "(log) Total Expenditure (Millions of US\\$)"),
  make_description_row(dem_cfscore,
    "Democrat CF Score"),
  make_description_row(rep_cfscore,
    "Republican CF Score"),
  make_description_row(ideological_distance,
    "Relative Republican Extremism"),
  make_description_row(adj_dem_presvote_advantage,
    "Adjusted Democratic Presidential Vote Advantage"),
  make_description_row(voting_eligible_population / 1e6,
    "Voting Eligible Population (in millions)"),
  make_description_row(open_seat,
    "Open Seat"),
  make_description_row(dem_incumbent,
    "Democratic Incumbent"),
  make_description_row(jungle,
    "Jungle Primary"),
  make_description_row(weird_race,
    "Top Two General Election/Other Candidate/Etc."),
  make_description_row(is.na(dem_vote + rep_vote), "Unopposed"))]
senate_data[special == 0, .N]
nrow(sen_X)
sink()

# Appendix E ----
#     Senate SVM ----
if (!file.exists(paste0(results_path, "senate_svm.RData"))) {
  senate_svm_model <- svm(x = sen_X, y = sen_y)
  cv_svm_sen <- cv_svm(sen_folds, sen_X, sen_y)
  senate_svm_data <- data.table(
    dem_spend_adv = sen_X[, "dem_spend_adv"],
    d_dem_spend_adv = ((
      predict(senate_svm_model, Xp1S) -
        predict(senate_svm_model, Xm1S)) /
        (2 * dxS)),
    incl = sen_X[, "bottom_tail"] + sen_X[, "top_tail"] == 0)
  save(senate_svm_model, cv_svm_sen, senate_svm_data,
    file = paste0(results_path, "senate_svm.RData"))
} else {
  load(paste0(results_path, "senate_svm.RData"))
}

#     Senate LASSO ----
if (!file.exists(paste0(results_path, "senate_lasso.RData"))) {
  senate_lasso_model <- glmnet(x = sen_X, y = sen_y,
    lambda = cv.glmnet(x = sen_X, y = sen_y)$lambda.1se)
  cv_lasso_sen <- cv_lasso(sen_folds, sen_X, sen_y)
  senate_lasso_data <- data.table(
    dem_spend_adv = sen_X[, "dem_spend_adv"],
    d_dem_spend_adv = ((
      predict(senate_lasso_model, Xp1S) -
        predict(senate_lasso_model, Xm1S)) /
        (2 * dxS))[, 1],
    incl = sen_X[, "bottom_tail"] + sen_X[, "top_tail"] == 0)
  save(senate_lasso_model, cv_lasso_sen, senate_lasso_data,
    file = paste0(results_path, "senate_lasso.RData"))
} else {
  load(paste0(results_path, "senate_lasso.RData"))
}

#     Senate Random Forests ----
if (!file.exists(paste0(results_path, "senate_rf.RData"))) {
  set.seed(520391846)
  senate_rf_model <- randomForest(x = sen_X_rf, y = sen_y[, 1], ntree = 1e5)
  senate_rf_data <- data.table(
    dem_spend_adv = sen_X[, "dem_spend_adv"],
    d_dem_spend_adv = ((
      predict(senate_rf_model, Xp1S) -
        predict(senate_rf_model, Xm1S)) /
        (2 * dxS)),
    incl = sen_X[, "bottom_tail"] + sen_X[, "top_tail"] == 0)
  cv_rf_sen <- cv_rf(sen_folds, sen_X_rf, sen_y, ntree = 1e5)
  save(senate_rf_model, cv_rf_sen, senate_rf_data,
    file = paste0(results_path, "senate_rf.RData"))
} else {
  load(paste0(results_path, "senate_rf.RData"))
}

#     House SVM ----
if (!file.exists(paste0(results_path, "house_svm.RData"))) {
  house_svm_model <- svm(x = hou_X, y = hou_y)
  cv_svm_hou <- cv_svm(hou_folds, hou_X, hou_y)
  house_svm_data <- data.table(
    dem_spend_adv = hou_X[, "dem_spend_adv"],
    d_dem_spend_adv = ((
      predict(house_svm_model, Xp1H) -
        predict(house_svm_model, Xm1H)) /
        (2 * dxH)),
    incl = hou_X[, "bottom_tail"] + hou_X[, "top_tail"] == 0)
  save(house_svm_model, cv_svm_hou, house_svm_data,
    file = paste0(results_path, "house_svm.RData"))
} else {
  load(paste0(results_path, "house_svm.RData"))
}

#     House LASSO ----
if (!file.exists(paste0(results_path, "house_lasso.RData"))) {
  house_lasso_model <- glmnet(x = hou_X, y = hou_y,
    lambda = cv.glmnet(x = hou_X, y = hou_y)$lambda.1se)
  cv_lasso_hou <- cv_lasso(hou_folds, hou_X, hou_y)
  house_lasso_data <- data.table(
    dem_spend_adv = hou_X[, "dem_spend_adv"],
    d_dem_spend_adv = ((
      predict(house_lasso_model, Xp1H) -
        predict(house_lasso_model, Xm1H)) /
        (2 * dxH))[, 1],
    incl = hou_X[, "bottom_tail"] + hou_X[, "top_tail"] == 0)
  save(house_lasso_model, cv_lasso_hou, house_lasso_data,
    file = paste0(results_path, "house_lasso.RData"))
} else {
  load(paste0(results_path, "house_lasso.RData"))
}

#     House Random Forests ----
if (!file.exists(paste0(results_path, "house_rf.RData"))) {
  set.seed(951358422)
  house_rf_model <- randomForest(x = hou_X_rf, y = hou_y[, 1], ntree = 1000)
  house_rf_data <- data.table(
    dem_spend_adv = hou_X[, "dem_spend_adv"],
    d_dem_spend_adv = ((
      predict(house_rf_model, Xp1H) -
        predict(house_rf_model, Xm1H)) /
        (2 * dxH)),
    incl = hou_X[, "bottom_tail"] + hou_X[, "top_tail"] == 0)
  cv_rf_hou <- cv_rf(hou_folds, hou_X, hou_y, ntree = 1000)
  save(house_rf_model, cv_rf_hou, house_rf_data,
    file = paste0(results_path, "house_rf.RData"))
} else {
  load(paste0(results_path, "house_rf.RData"))
}

# ranges of random forest marginal effect estimates ----
# these should be bounded by [-1,1], but they are not
house_rf_data[, range(d_dem_spend_adv)]
senate_rf_data[, range(d_dem_spend_adv)]
# > 2% of ME estimates are out of bounds in the House
house_rf_data[, mean(abs(d_dem_spend_adv) > 1)]
senate_rf_data[, mean(abs(d_dem_spend_adv) > 1)]

# this leads to astronomical estimates for the effectiveness of spending
ggplot(house_rf_data[dem_spend_adv != 0],
  aes(dem_spend_adv, d_dem_spend_adv)) + geom_smooth() +
  theme(text = element_text(family = "sans"))
ggplot(senate_rf_data[dem_spend_adv != 0],
  aes(dem_spend_adv, d_dem_spend_adv)) + geom_smooth() +
  theme(text = element_text(family = "sans"))

# when we drop those outliers, we get estimates that match other methods
ggplot(house_rf_data[abs(d_dem_spend_adv) < 1 & dem_spend_adv != 0],
  aes(dem_spend_adv, d_dem_spend_adv)) + geom_smooth() +
  theme(text = element_text(family = "sans"))
ggplot(senate_rf_data[abs(d_dem_spend_adv) < 1 & dem_spend_adv != 0],
  aes(dem_spend_adv, d_dem_spend_adv)) + geom_smooth() +
  theme(text = element_text(family = "sans"))

#     Make CV/OOS Table (S5) ----
oos_results <- CJ(
  chamber = c("House", "Senate"),
  method = c("KRLS", "SVM", "LASSO", "RF"),
  measure = c("RMSE", "R2"),
  value = NA_real_)
oos_results[chamber == "Senate" & method == "KRLS" & measure == "RMSE",
  value := mean(cv_rmse(senate_krls_cv, sen_folds, sen_y))]
oos_results[chamber == "Senate" & method == "KRLS" & measure == "R2",
  value := mean(cv_R2(senate_krls_cv, sen_folds, sen_y))]
oos_results[chamber == "Senate" & method == "RF" & measure == "RMSE",
  value := mean(cv_rmse(cv_rf_sen, sen_folds, sen_y))]
oos_results[chamber == "Senate" & method == "RF" & measure == "R2",
  value := mean(cv_R2(cv_rf_sen, sen_folds, sen_y))]
oos_results[chamber == "Senate" & method == "SVM" & measure == "RMSE",
  value := mean(cv_rmse(cv_svm_sen, sen_folds, sen_y))]
oos_results[chamber == "Senate" & method == "SVM" & measure == "R2",
  value := mean(cv_R2(cv_svm_sen, sen_folds, sen_y))]
oos_results[chamber == "Senate" & method == "LASSO" & measure == "RMSE",
  value := mean(cv_rmse(cv_lasso_sen, sen_folds, sen_y))]
oos_results[chamber == "Senate" & method == "LASSO" & measure == "R2",
  value := mean(cv_R2(cv_lasso_sen, sen_folds, sen_y))]
oos_results[chamber == "House" & method == "KRLS" & measure == "RMSE",
  value := mean(cv_rmse(house_krls_cv, hou_folds, hou_y))]
oos_results[chamber == "House" & method == "KRLS" & measure == "R2",
  value := mean(cv_R2(house_krls_cv, hou_folds, hou_y))]
oos_results[chamber == "House" & method == "RF" & measure == "RMSE",
  value := mean(cv_rmse(cv_rf_hou, hou_folds, hou_y))]
oos_results[chamber == "House" & method == "RF" & measure == "R2",
  value := mean(cv_R2(cv_rf_hou, hou_folds, hou_y))]
oos_results[chamber == "House" & method == "SVM" & measure == "RMSE",
  value := mean(cv_rmse(cv_svm_hou, hou_folds, hou_y))]
oos_results[chamber == "House" & method == "SVM" & measure == "R2",
  value := mean(cv_R2(cv_svm_hou, hou_folds, hou_y))]
oos_results[chamber == "House" & method == "LASSO" & measure == "RMSE",
  value := mean(cv_rmse(cv_lasso_hou, hou_folds, hou_y))]
oos_results[chamber == "House" & method == "LASSO" & measure == "R2",
  value := mean(cv_R2(cv_lasso_hou, hou_folds, hou_y))]
oos_results[, value := round(value, 3)]
sink("tables/table-s5-top.txt")
dcast(oos_results[measure == "RMSE"], method ~ chamber, value.vars = "value")[
  c(1, 3, 4, 2)]
sink()
sink("tables/table-s5-bottom.txt")
dcast(oos_results[measure == "R2"], method ~ chamber, value.vars = "value")[
  c(1, 3, 4, 2)]
sink()

#     Make average marginal effects Table (S6) ----
ame_robustness <- rbind(
  senate_krls_results[bottom_tail + top_tail == 0, .(
    method = "KRLS", chamber = "Senate",
    Estimate = mean(d_dem_spend_adv), SE = calc_sem(d_dem_spend_adv, .N),
    pval = calc_p(d_dem_spend_adv, .N))],
  senate_rf_data[incl == TRUE, .(
    method = "RF", chamber = "Senate",
    Estimate = mean(d_dem_spend_adv), SE = calc_sem(d_dem_spend_adv, .N),
    pval = calc_p(d_dem_spend_adv, .N))],
  senate_svm_data[incl == TRUE, .(
    method = "SVM", chamber = "Senate",
    Estimate = mean(d_dem_spend_adv), SE = calc_sem(d_dem_spend_adv, .N),
    pval = calc_p(d_dem_spend_adv, .N))],
  senate_lasso_data[incl == TRUE, .(
    method = "LASSO", chamber = "Senate",
    Estimate = mean(d_dem_spend_adv), SE = calc_sem(d_dem_spend_adv, .N),
    pval = calc_p(d_dem_spend_adv, .N))],
  house_krls_results[bottom_tail + top_tail == 0, .(
    method = "KRLS", chamber = "House",
    Estimate = mean(d_dem_spend_adv), SE = calc_sem(d_dem_spend_adv, .N),
    pval = calc_p(d_dem_spend_adv, .N))],
  house_rf_data[incl == TRUE, .(
    method = "RF", chamber = "House",
    Estimate = mean(d_dem_spend_adv), SE = calc_sem(d_dem_spend_adv, .N),
    pval = calc_p(d_dem_spend_adv, .N))],
  house_svm_data[incl == TRUE, .(
    method = "SVM", chamber = "House",
    Estimate = mean(d_dem_spend_adv), SE = calc_sem(d_dem_spend_adv, .N),
    pval = calc_p(d_dem_spend_adv, .N))],
  house_lasso_data[incl == TRUE, .(
    method = "LASSO", chamber = "House",
    Estimate = mean(d_dem_spend_adv), SE = calc_sem(d_dem_spend_adv, .N),
    pval = calc_p(d_dem_spend_adv, .N))]
)
ames_robustness_summ <- ame_robustness[, .(method, chamber,
  Estimate = round(Estimate, 3), SE = round(SE, 3),
  pval = round(pval, 3))]
ames_robustness_summ[chamber == "House"]
ames_robustness_summ[chamber == "Senate"]
#     LASSO with polynomial in spending advantage ----
if (!file.exists(paste0(results_path, "sen_lasso_poly.RData"))) {
  sen_X_poly <- cbind(sen_X,
    sen_X[, "dem_spend_adv"]^2, sen_X[, "dem_spend_adv"]^3,
    sen_X[, "dem_spend_adv"]^4,
    sen_X[, "bottom_tail_DSA"]^2, sen_X[, "bottom_tail_DSA"]^3,
    sen_X[, "bottom_tail_DSA"]^4,
    sen_X[, "top_tail_DSA"]^2, sen_X[, "top_tail_DSA"]^3,
    sen_X[, "top_tail_DSA"]^4)
  senate_lasso_poly_model <- glmnet(x = sen_X_poly, y = sen_y,
    lambda = cv.glmnet(x = sen_X_poly, y = sen_y)$lambda.1se)
  cv_lasso_poly_sen <- cv_lasso(sen_folds, sen_X_poly, sen_y)
  Xp1S_poly <- cbind(Xp1S,
    Xp1S[, "dem_spend_adv"]^2, Xp1S[, "dem_spend_adv"]^3,
    Xp1S[, "dem_spend_adv"]^4,
    Xp1S[, "bottom_tail_DSA"]^2, Xp1S[, "bottom_tail_DSA"]^3,
    Xp1S[, "bottom_tail_DSA"]^4,
    Xp1S[, "top_tail_DSA"]^2, Xp1S[, "top_tail_DSA"]^3,
    Xp1S[, "top_tail_DSA"]^4)
  Xm1S_poly <- cbind(Xm1S,
    Xm1S[, "dem_spend_adv"]^2, Xm1S[, "dem_spend_adv"]^3,
    Xm1S[, "dem_spend_adv"]^4,
    Xm1S[, "bottom_tail_DSA"]^2, Xm1S[, "bottom_tail_DSA"]^3,
    Xm1S[, "bottom_tail_DSA"]^4,
    Xm1S[, "top_tail_DSA"]^2, Xm1S[, "top_tail_DSA"]^3,
    Xm1S[, "top_tail_DSA"]^4)
  senate_lasso_poly_data <- data.table(
    dem_spend_adv = sen_X_poly[, "dem_spend_adv"],
    d_dem_spend_adv = ((
      predict(senate_lasso_poly_model, Xp1S_poly) -
        predict(senate_lasso_poly_model, Xm1S_poly)) /
        (2 * dxS))[, 1],
    incl = sen_X_poly[, "bottom_tail"] + sen_X_poly[, "top_tail"] == 0)
  save(senate_lasso_poly_model, cv_lasso_poly_sen, senate_lasso_poly_data,
    file = paste0(results_path, "sen_lasso_poly.RData"))
} else {
  load(paste0(results_path, "sen_lasso_poly.RData"))
}
if (!file.exists(paste0(results_path, "house_lasso_poly.RData"))) {
  hou_X_poly <- cbind(hou_X,
    hou_X[, "dem_spend_adv"]^2, hou_X[, "dem_spend_adv"]^3,
    hou_X[, "dem_spend_adv"]^4,
    hou_X[, "bottom_tail_DSA"]^2, hou_X[, "bottom_tail_DSA"]^3,
    hou_X[, "bottom_tail_DSA"]^4,
    hou_X[, "top_tail_DSA"]^2, hou_X[, "top_tail_DSA"]^3,
    hou_X[, "top_tail_DSA"]^4)
  house_lasso_poly_model <- glmnet(x = hou_X_poly, y = hou_y,
    lambda = cv.glmnet(x = hou_X_poly, y = hou_y)$lambda.1se)
  cv_lasso_poly_hou <- cv_lasso(hou_folds, hou_X_poly, hou_y)
  Xp1H_poly <- cbind(Xp1H,
    Xp1H[, "dem_spend_adv"]^2, Xp1H[, "dem_spend_adv"]^3,
    Xp1H[, "dem_spend_adv"]^4,
    Xp1H[, "bottom_tail_DSA"]^2, Xp1H[, "bottom_tail_DSA"]^3,
    Xp1H[, "bottom_tail_DSA"]^4,
    Xp1H[, "top_tail_DSA"]^2, Xp1H[, "top_tail_DSA"]^3,
    Xp1H[, "top_tail_DSA"]^4)
  Xm1H_poly <- cbind(Xm1H,
    Xm1H[, "dem_spend_adv"]^2, Xm1H[, "dem_spend_adv"]^3,
    Xm1H[, "dem_spend_adv"]^4,
    Xm1H[, "bottom_tail_DSA"]^2, Xm1H[, "bottom_tail_DSA"]^3,
    Xm1H[, "bottom_tail_DSA"]^4,
    Xm1H[, "top_tail_DSA"]^2, Xm1H[, "top_tail_DSA"]^3,
    Xm1H[, "top_tail_DSA"]^4)
  house_lasso_poly_data <- data.table(
    dem_spend_adv = hou_X_poly[, "dem_spend_adv"],
    d_dem_spend_adv = ((
      predict(house_lasso_poly_model, Xp1H_poly) -
        predict(house_lasso_poly_model, Xm1H_poly)) /
        (2 * dxH))[, 1],
    incl = hou_X_poly[, "bottom_tail"] + hou_X_poly[, "top_tail"] == 0)
  save(house_lasso_poly_model, cv_lasso_poly_hou, house_lasso_poly_data,
    file = paste0(results_path, "house_lasso_poly.RData"))
} else {
  load(paste0(results_path, "house_lasso_poly.RData"))
}
round(mean(cv_rmse(cv_lasso_poly_sen, sen_folds, sen_y)), 3)
round(mean(cv_R2(cv_lasso_poly_sen, sen_folds, sen_y)), 2)

round(mean(cv_rmse(cv_lasso_poly_hou, hou_folds, hou_y)), 3)
round(mean(cv_R2(cv_lasso_poly_hou, hou_folds, hou_y)), 2)


#     Make Figure A1 (robustness of Figure 1) ----
oos_data <- rbind(
  house_rf_data[incl == TRUE,
    .(dem_spend_adv, d_dem_spend_adv, Method = "RF", chamber = "House")],
  house_svm_data[incl == TRUE,
    .(dem_spend_adv, d_dem_spend_adv, Method = "SVM", chamber = "House")],
  house_lasso_data[incl == TRUE,
    .(dem_spend_adv, d_dem_spend_adv = round(d_dem_spend_adv, 4),
      Method = "LASSO", chamber = "House")],
  house_krls_results[bottom_tail+top_tail == 0,
    .(dem_spend_adv, d_dem_spend_adv, Method = "KRLS", chamber = "House")],
  senate_rf_data[incl == TRUE,
    .(dem_spend_adv, d_dem_spend_adv, Method = "RF", chamber = "Senate")],
  senate_svm_data[incl == TRUE,
    .(dem_spend_adv, d_dem_spend_adv, Method = "SVM", chamber = "Senate")],
  senate_lasso_data[incl == TRUE,
    .(dem_spend_adv,  d_dem_spend_adv = round(d_dem_spend_adv, 4),
      Method = "LASSO", chamber = "Senate")],
  senate_krls_results[bottom_tail+top_tail == 0,
    .(dem_spend_adv, d_dem_spend_adv, Method = "KRLS", chamber = "Senate")]
)
fig_oos <- ggplot(oos_data[Method != "LASSO"], aes(dem_spend_adv, d_dem_spend_adv,
  color = Method, linetype = Method)) +
  geom_smooth(se = FALSE, method = "loess", span = 1) +
  facet_wrap(~chamber, scales = "free") +
  xlab("Democratic Spending Advantage") +
  theme(legend.position = "bottom") +
  ylab("Estimated Average Marginal Effects") +
  scale_linetype_manual(values = c(
    "KRLS" = "solid", "SVM" = "dotted", "RF" = "longdash"
  ))

cairo_pdf(paste0("figures/fig-oos-", Sys.Date(), ".pdf"),
  width = 6.5, height = 4)
h <- .15
h2 <- .03
suppressWarnings(
  cowplot::ggdraw() +
    cowplot::draw_plot(fig_oos + theme(
      plot.title = element_text(size = 12),
      text = element_text(family = "CM Roman")),
      x = 0, y = .35 * h,  height = 1 - 1.5 * h, width = 1) +
    cowplot::draw_label(
      paste0("Robustness of Nonlinear Effect Estimates"),
      size = 18, x = .5, y = 1 - .8 * h, fontface = "bold",
      fontfamily = "CM Roman")
)
dev.off()

# Appendix F ----
#     House KRLS w/ lagged outcome ----
if (!file.exists(paste0(results_path, "house_krls_lag.RData"))) {
  house_krls_model_lag <- bigKRLS::bigKRLS(X = hou_X_lag,
    y = hou_y_lag, noisy = TRUE)
  house_krls_results_lag <- as.data.table(cbind(hou_X_lag,
    house_krls_model_lag$derivatives))
  setnames(house_krls_results_lag, c(colnames(hou_X_lag),
    paste0("d_", colnames(hou_X_lag))))
  house_krls_sum_lag <- summary(house_krls_model_lag)
  house_krls_sum_lag$ttests <- as.data.table(house_krls_sum_lag$ttests)
  house_krls_sum_lag$percentiles <-
    as.data.table(house_krls_sum_lag$percentiles)
  save.bigKRLS(house_krls_model_lag,
    model_subfolder_name = paste0(results_path, "house_krls_model_lag"))
  save(house_krls_results_lag, house_krls_sum_lag,
    file = paste0(results_path, "house_krls_lag.RData"))
} else {
  load.bigKRLS(paste0(results_path, "house_krls_model_lag"),
    newname = "house_krls_model_lag")
  load(paste0(results_path, "house_krls_lag.RData"))
}

hou_index_lag <- rbind(
  c(1, "Democratic Spending Advantage"),
  c(2, "$\\text{log}$(Total Spending)"),
  c(3, "Democratic Presidential Vote Advantage"),
  c(4, "Relative Republican Extremism"),
  c(5, "Democrat Candidate's CF Score"),
  c(6, "Republican Candidate's CF Score"),
  c(7, "Dem. Inc./Low Qual. Chall."),
  c(8, "Dem. Inc./High Qual. Chall."),
  c(9, "Rep Inc./Low Qual. Chall."),
  c(10, "Rep Inc./High Qual. Chall."),
  c(11, "Year = 1984"),
  c(12, "Year = 1986"),
  c(13, "Year = 1988"),
  c(14, "Year = 1990"),
  c(15, "Year = 1994"),
  c(16, "Year = 1996"),
  c(17, "Year = 1998"),
  c(18, "Year = 2000"),
  c(19, "Year = 2004"),
  c(20, "Year = 2006"),
  c(21, "Year = 2008"),
  c(22, "Year = 2010"),
  c(23, "Year = 2014"),
  c(24, "Year = 2016"),
  c(25, "Year = 2018"),
  c(26, "Bottom 5\\%"),
  c(27, "Top 5\\%"),
  c(30, "Middle 90\\%"),
  c(28, "Bottom 5\\% $\\times$ Dem. Spending Adv."),
  c(29, "Top 5\\% $\\times$ Dem. Spending Adv."),
  c(31, "Lagged Outcome")
)
for (i in 1:nrow(hou_index_lag)) {
  make_outcome_row(house_krls_sum_lag$ttests[as.numeric(hou_index_lag[i, 1]), ],
    hou_index_lag[i, 2])
}

#     Senate KRLS w/ lagged outcome ----
if (!file.exists(paste0(results_path, "senate_krls_model_lag"))) {
  senate_krls_model_lag <- bigKRLS::bigKRLS(X = sen_X_lag,
    y = sen_y_lag, noisy = TRUE)
  save.bigKRLS(senate_krls_model_lag,
    model_subfolder_name = paste0(results_path, "senate_krls_model_lag"))
} else {
  load.bigKRLS(paste0(results_path, "senate_krls_model_lag"),
    newname = "senate_krls_model_lag")
}
senate_krls_sum_lag <- summary(senate_krls_model_lag)
senate_krls_sum_lag$ttests <- as.data.table(senate_krls_sum_lag$ttests)
senate_krls_sum_lag$percentiles <-
  as.data.table(senate_krls_sum_lag$percentiles)
sen_index_lag <- rbind(
  c(1, "Democratic Spending Advantage"),
  c(2, "$\\text{log}$(Total Spending)"),
  c(3, "Democratic Presidential Vote Advantage"),
  c(4, "Adjusted Democratic Presidential Vote Advantage"),
  c(5, "$\\text{log}$($n$ Votes for Democratic Presidential Candidate)"),
  c(6, "$\\text{log}$($n$ Votes for Republican Presidential Candidate)"),
  c(7, "Relative Republican Extremism"),
  c(8, "Democrat Candidate's CF Score"),
  c(9, "Republican Candidate's CF Score"),
  c(10, "Democratic Incumbent"),
  c(11, "$\\text{log}$(Voting Eligible Population)"),
  c(12, "Year = 1986"),
  c(13, "Year = 1988"),
  c(14, "Year = 1990"),
  c(15, "Year = 1992"),
  c(16, "Year = 1994"),
  c(17, "Year = 1996"),
  c(18, "Year = 1998"),
  c(19, "Year = 2000"),
  c(20, "Year = 2002"),
  c(21, "Year = 2004"),
  c(22, "Year = 2006"),
  c(23, "Year = 2008"),
  c(24, "Year = 2010"),
  c(25, "Year = 2012"),
  c(26, "Year = 2014"),
  c(27, "Year = 2016"),
  c(28, "Year = 2018"),
  c(29, "Bottom 5\\%"),
  c(30, "Top 5\\%"),
  c(33, "Middle 90\\%"),
  c(31, "Bottom 5\\% $\\times$ Dem. Spending Adv."),
  c(32, "Top 5\\% $\\times$ Dem. Spending Adv."),
  c(34, "Lagged Outcome")
)
for (i in 1:nrow(sen_index_lag)) {
  make_outcome_row(senate_krls_sum_lag$ttests[
    as.numeric(sen_index_lag[i, 1]), ], sen_index_lag[i, 2])
}

#     Senate and House KRLS w/ lagged outcome + lagged spending ----
if (!file.exists(paste0(results_path, "senate_krls_model_lag2"))) {
  sen_X_lag2 <- cbind(sen_X_lag, data.matrix(senate_data_lag[, lag_DSA]))
  colnames(sen_X_lag2)[ncol(sen_X_lag2)] <- "lag_DSA"
  ok <- complete.cases(sen_X_lag2)
  sen_X_lag2 <- sen_X_lag2[ok, ]
  sen_y_lag2 <- sen_y_lag[ok, , drop = FALSE]
  senate_krls_model_lag2 <- bigKRLS::bigKRLS(X = sen_X_lag2,
    y = sen_y_lag2, noisy = TRUE)
  summary(senate_krls_model_lag2)

  hou_X_lag2 <- cbind(hou_X_lag, data.matrix(house_data_lag[, lag_DSA]))
  colnames(hou_X_lag2)[ncol(hou_X_lag2)] <- "lag_DSA"
  ok <- complete.cases(hou_X_lag2)
  hou_X_lag2 <- hou_X_lag2[ok, ]
  hou_y_lag2 <- hou_y_lag[ok, , drop = FALSE]
  house_krls_model_lag2 <- bigKRLS::bigKRLS(X = hou_X_lag2,
    y = hou_y_lag2, noisy = TRUE)
  summary(house_krls_model_lag2)

  save.bigKRLS(house_krls_model_lag2,
    model_subfolder_name = paste0(results_path, "house_krls_model_lag2"))
  save.bigKRLS(senate_krls_model_lag2,
    model_subfolder_name = paste0(results_path, "senate_krls_model_lag2"))
} else {
  load.bigKRLS(paste0(results_path, "house_krls_model_lag2"),
    newname = "house_krls_model_lag2")
  load.bigKRLS(paste0(results_path, "senate_krls_model_lag2"),
    newname = "senate_krls_model_lag2")
}
# Appendix G ----
#     House KRLS w/ dichotomous outcomes ----
if (!file.exists(paste0(results_path, "house_krls_model_01"))) {
  house_krls_model_01 <- bigKRLS::bigKRLS(X = hou_X, y = hou_y_01, noisy = TRUE)
  save.bigKRLS(house_krls_model_01,
    model_subfolder_name = paste0(results_path, "house_krls_model_01"))
  house_krls_results_01 <- as.data.table(cbind(hou_X,
    house_krls_model_01$derivatives))
  setnames(house_krls_results_01, c(colnames(hou_X),
    paste0("d_", colnames(hou_X))))
  house_krls_sum_01 <- summary(house_krls_model_01)
  house_krls_sum_01$ttests <- as.data.table(house_krls_sum_01$ttests)
  house_krls_sum_01$percentiles <- as.data.table(house_krls_sum_01$percentiles)
  save(house_krls_results_01, house_krls_sum_01,
    file = paste0(results_path, "house_krls_01.RData"))
} else {
  load.bigKRLS(paste0(results_path, "house_krls_model_01"))
  load(paste0(results_path, "house_krls_01.RData"))
}
for (i in 1:nrow(hou_index)) {
  make_outcome_row(house_krls_sum_01$ttests[
    as.numeric(hou_index[i, 1]), ], hou_index[i, 2])
}

#     Senate KRLS w/ dichotomous outcomes ----
if (!file.exists(paste0(results_path, "senate_krls_model_01"))) {
  senate_krls_model_01 <- bigKRLS::bigKRLS(X = sen_X, y = sen_y_01,
    noisy = TRUE)
  save.bigKRLS(senate_krls_model_01,
    model_subfolder_name = paste0(results_path, "senate_krls_model_01"))
} else {
  load.bigKRLS(paste0(results_path, "senate_krls_model_01"),
    newname = "senate_krls_model_01")
}
senate_krls_results_01 <- as.data.table(cbind(sen_X,
  senate_krls_model_01$derivatives))
setnames(senate_krls_results_01, c(colnames(sen_X),
  paste0("d_", colnames(sen_X))))
senate_krls_sum_01 <- summary(senate_krls_model_01)
senate_krls_sum_01$ttests <- as.data.table(senate_krls_sum_01$ttests)
senate_krls_sum_01$percentiles <-
  as.data.table(senate_krls_sum_01$percentiles)
# make table of results for dichot model
for (i in 1:nrow(sen_index)) {
  make_outcome_row(senate_krls_sum_01$ttests[
    as.numeric(sen_index[i, 1]), ], sen_index[i, 2])
}


# Appendix H ----
#   SVM results ----
house_svm_model <- svm(x = hou_X, y = hou_y)
house_svm_data <- data.table(
  dem_spend_adv = hou_X[, "dem_spend_adv"],
  d_dem_spend_adv = ((
    predict(house_svm_model, Xp1H) -
      predict(house_svm_model, Xm1H)) /
      (2 * dxH)),
  incl = hou_X[, "bottom_tail"] + hou_X[, "top_tail"] == 0)
senate_svm_model <- svm(x = sen_X, y = sen_y)
senate_svm_data <- data.table(
  dem_spend_adv = sen_X[, "dem_spend_adv"],
  d_dem_spend_adv = ((
    predict(senate_svm_model, Xp1S) -
      predict(senate_svm_model, Xm1S)) /
      (2 * dxS)),
  incl = sen_X[, "bottom_tail"] + sen_X[, "top_tail"] == 0)

loess_house <- house_svm_data[incl == TRUE,
  loess(d_dem_spend_adv ~ dem_spend_adv)]
integrate(function(x) predict(loess_house, x),
  lower = house_svm_data[incl == TRUE, min(dem_spend_adv)],
  upper = house_svm_data[incl == TRUE, max(dem_spend_adv)])

loess_senate <- senate_svm_data[incl == TRUE,
  loess(d_dem_spend_adv ~ dem_spend_adv)]
integrate(function(x) predict(loess_senate, x),
  lower = senate_svm_data[incl == TRUE, min(dem_spend_adv)],
  upper = senate_svm_data[incl == TRUE, max(dem_spend_adv)])
#     Senate SVM counterfactuals ----
# estimate models to nonparametric bootstrapped datasets
if (!file.exists(paste0(results_path, "senate_boot_svm_models.RData"))) {
  set.seed(990264419)
  senate_boot_svm_models <- lapply(1:1000, function(i) {
    ok <- sample(nrow(sen_X), nrow(sen_X), replace = TRUE)
    svm(x = sen_X[ok, ], y = sen_y[ok])
  })
  save(senate_boot_svm_models,
    file = paste0(results_path, "senate_boot_svm_models.RData"))
} else {
  load(paste0(results_path, "senate_boot_svm_models.RData"))
}
# calc marginal effects on bootstrapped datasets,
#   then estimate loess curves of those marginal effects
if (!file.exists(paste0(results_path, "senate_boot_svm_cond_ame.RData"))) {
  set.seed(990264419)
  senate_boot_svm_cond_ame <- do.call(cbind, lapply(1:1000, function(i) {
    ok <- sample(nrow(sen_X), nrow(sen_X), replace = TRUE)
    xS_boot <- sen_X[ok, 1]
    dxS_boot <- setstep(xS_boot)
    Xp1S_boot <- make_plus_dx_dataset(sen_X[ok, ], dxS_boot, "sen")
    Xm1S_boot <- make_plus_dx_dataset(sen_X[ok, ], -dxS_boot, "sen")
    senate_boot_svm_data <- data.table(
      dem_spend_adv = sen_X[ok, "dem_spend_adv"],
      d_dem_spend_adv = ((
        predict(senate_boot_svm_models[[i]], Xp1S_boot) -
          predict(senate_boot_svm_models[[i]], Xm1S_boot)) /
          (2 * dxS_boot)),
      incl = sen_X[ok, "bottom_tail"] + sen_X[ok, "top_tail"] == 0)
    predict(
      senate_boot_svm_data[incl == TRUE, loess(d_dem_spend_adv ~ dem_spend_adv)],
      sen_X[sen_X[, "bottom_tail"] + sen_X[, "top_tail"] == 0, ])
  }))
  save(senate_boot_svm_cond_ame,
    file = paste0(results_path, "senate_boot_svm_cond_ame.RData"))
} else {
  load(paste0(results_path, "senate_boot_svm_cond_ame.RData"))
}
# estimate counterfactuals
if (!file.exists(paste0(results_path, "senate_svm_counterfactuals.RData"))) {
  n_boots <- 1000
  n_obs <- nrow(sen_X)
  possible_values <- c(NA, 0,
    sen_counterfactual_vals[2], sen_counterfactual_vals[1])
  senate_svm_counterfactuals <- CJ(boot = 1:n_boots,
    dem_spend_adv = possible_values, index = senate_data[, .I])
  senate_svm_counterfactuals[, actual := as.numeric(is.na(dem_spend_adv))]
  senate_svm_counterfactuals[is.na(dem_spend_adv),
    dem_spend_adv := sen_X[index, 1]]
  year <-
    1980 * sen_X[, "y80"] + 1982 * sen_X[, "y82"] + 1984 * sen_X[, "y84"] +
    1986 * sen_X[, "y86"] + 1988 * sen_X[, "y88"] + 1990 * sen_X[, "y90"] +
    1992 * sen_X[, "y92"] + 1994 * sen_X[, "y94"] + 1996 * sen_X[, "y96"] +
    1998 * sen_X[, "y98"] + 2000 * sen_X[, "y00"] + 2002 * sen_X[, "y02"] +
    2004 * sen_X[, "y04"] + 2006 * sen_X[, "y06"] + 2008 * sen_X[, "y08"] +
    2010 * sen_X[, "y10"] + 2012 * sen_X[, "y12"] + 2014 * sen_X[, "y14"] +
    2016 * sen_X[, "y16"] + 2018 * sen_X[, "y18"]
  senate_svm_counterfactuals[, year := year[index]]
  # actual
  sen_pred_X <- sen_X
  pred <- do.call(rbind, lapply(senate_boot_svm_models, predict,
    newdata = sen_pred_X))
  senate_svm_counterfactuals[actual == 1,
    boot_y := pred[cbind(boot,index)]]
  # zero spending adv
  sen_pred_X <- sen_X
  sen_pred_X[, 1] <- 0
  sen_pred_X[, 2] <- min(sen_X[, "log_total_spending"])
  pred <- do.call(rbind, lapply(senate_boot_svm_models, predict,
    newdata = sen_pred_X))
  senate_svm_counterfactuals[actual == 0 & dem_spend_adv == 0,
    boot_y := pred[cbind(boot,index)]]
  # max dem spending
  sen_pred_X <- sen_X
  sen_pred_X[, 1] <- sen_counterfactual_vals[2]
  sen_pred_X[, 2] <- senate_data[, log10(total_spending_dem_max)]
  pred <- do.call(rbind, lapply(senate_boot_svm_models, predict,
    newdata = sen_pred_X))
  senate_svm_counterfactuals[actual == 0 &
      dem_spend_adv == sen_counterfactual_vals[2],
    boot_y := pred[cbind(boot,index)]]
  # max rep spending
  sen_pred_X <- sen_X
  sen_pred_X[, 1] <- sen_counterfactual_vals[1]
  sen_pred_X[, 2] <- senate_data[, log10(total_spending_rep_max)]
  pred <- do.call(rbind, lapply(senate_boot_svm_models, predict,
    newdata = sen_pred_X))
  senate_svm_counterfactuals[actual == 0 &
      dem_spend_adv == sen_counterfactual_vals[1],
    boot_y := pred[cbind(boot,index)]]
  save(senate_svm_counterfactuals,
    file = paste0(results_path, "senate_svm_counterfactuals.RData"))
} else {
  load(paste0(results_path, "senate_svm_counterfactuals.RData"))
}
#     House SVM counterfactuals ----
# estimate models to nonparametric bootstrapped datasets
if (!file.exists(paste0(results_path, "house_boot_svm_models.RData"))) {
  set.seed(990264419)
  house_boot_svm_models <- lapply(1:1000, function(i) {
    cat("\r", i)
    ok <- sample(nrow(hou_X), nrow(hou_X), replace = TRUE)
    svm(x = hou_X[ok, ], y = hou_y[ok])
  })
  save(house_boot_svm_models,
    file = paste0(results_path, "house_boot_svm_models.RData"))
} else {
  load(paste0(results_path, "house_boot_svm_models.RData"))
}
# calc marginal effects on bootstrapped datasets,
#   then estimate loess curves of those marginal effects
if (!file.exists(paste0(results_path, "house_boot_svm_cond_ame.RData"))) {
  set.seed(54266132)
  house_boot_svm_cond_ame <- do.call(cbind, lapply(1:1000, function(i) {
    cat("\r", i)
    ok <- sample(nrow(hou_X), nrow(hou_X), replace = TRUE)
    xH_boot <- hou_X[ok, 1]
    dxH_boot <- setstep(xH_boot)
    Xp1H_boot <- make_plus_dx_dataset(hou_X[ok, ], dxH_boot, "hou")
    Xm1H_boot <- make_plus_dx_dataset(hou_X[ok, ], -dxH_boot, "hou")
    house_boot_svm_data <- data.table(
      dem_spend_adv = hou_X[ok, "dem_spend_adv"],
      d_dem_spend_adv = ((
        predict(house_boot_svm_models[[i]], Xp1H_boot) -
          predict(house_boot_svm_models[[i]], Xm1H_boot)) /
          (2 * dxH_boot)),
      incl = hou_X[ok, "bottom_tail"] + hou_X[ok, "top_tail"] == 0)
    predict(
      house_boot_svm_data[incl == TRUE, loess(d_dem_spend_adv ~ dem_spend_adv)],
      hou_X[hou_X[, "bottom_tail"] + hou_X[, "top_tail"] == 0, ])
  }))
  save(house_boot_svm_cond_ame,
    file = paste0(results_path, "house_boot_svm_cond_ame.RData"))
} else {
  load(paste0(results_path, "house_boot_svm_cond_ame.RData"))
}
# estimate counterfactuals
if (!file.exists(paste0(results_path, "house_svm_counterfactuals.RData"))) {
  n_boots <- 1000
  n_obs <- nrow(hou_X)
  possible_values <- c(NA, 0,
    hou_counterfactual_vals[2], hou_counterfactual_vals[1])
  house_svm_counterfactuals <- CJ(boot = 1:n_boots,
    dem_spend_adv = possible_values, index = house_data[, .I])
  house_svm_counterfactuals[, actual := as.numeric(is.na(dem_spend_adv))]
  house_svm_counterfactuals[is.na(dem_spend_adv),
    dem_spend_adv := hou_X[index, 1]]
  year <-
    1980 * hou_X[, "y80"] + 1982 * hou_X[, "y82"] + 1984 * hou_X[, "y84"] +
    1986 * hou_X[, "y86"] + 1988 * hou_X[, "y88"] + 1990 * hou_X[, "y90"] +
    1992 * hou_X[, "y92"] + 1994 * hou_X[, "y94"] + 1996 * hou_X[, "y96"] +
    1998 * hou_X[, "y98"] + 2000 * hou_X[, "y00"] + 2002 * hou_X[, "y02"] +
    2004 * hou_X[, "y04"] + 2006 * hou_X[, "y06"] + 2008 * hou_X[, "y08"] +
    2010 * hou_X[, "y10"] + 2012 * hou_X[, "y12"] + 2014 * hou_X[, "y14"] +
    2016 * hou_X[, "y16"] + 2018 * hou_X[, "y18"]
  house_svm_counterfactuals[, year := year[index]]
  # actual
  hou_pred_X <- hou_X
  pred <- do.call(rbind, lapply(house_boot_svm_models, predict,
    newdata = hou_pred_X))
  house_svm_counterfactuals[actual == 1,
    boot_y := pred[cbind(boot,index)]]
  # zero spending adv
  hou_pred_X <- hou_X
  hou_pred_X[, 1] <- 0
  hou_pred_X[, 2] <- min(hou_X[, "log_total_spending"])
  pred <- do.call(rbind, lapply(house_boot_svm_models, predict,
    newdata = hou_pred_X))
  house_svm_counterfactuals[actual == 0 & dem_spend_adv == 0,
    boot_y := pred[cbind(boot,index)]]
  # max dem spending
  hou_pred_X <- hou_X
  hou_pred_X[, 1] <- hou_counterfactual_vals[2]
  hou_pred_X[, 2] <- house_data[, log10(total_spending_dem_max)]
  pred <- do.call(rbind, lapply(house_boot_svm_models, predict,
    newdata = hou_pred_X))
  house_svm_counterfactuals[actual == 0 &
      dem_spend_adv == hou_counterfactual_vals[2],
    boot_y := pred[cbind(boot,index)]]
  # max rep spending
  hou_pred_X <- hou_X
  hou_pred_X[, 1] <- hou_counterfactual_vals[1]
  hou_pred_X[, 2] <- house_data[, log10(total_spending_rep_max)]
  pred <- do.call(rbind, lapply(house_boot_svm_models, predict,
    newdata = hou_pred_X))
  house_svm_counterfactuals[actual == 0 &
      dem_spend_adv == hou_counterfactual_vals[1],
    boot_y := pred[cbind(boot,index)]]
  save(house_svm_counterfactuals,
    file = paste0(results_path, "house_svm_counterfactuals.RData"))
} else {
  load(paste0(results_path, "house_svm_counterfactuals.RData"))
}
#     Senate SVM figure data ----
if (!file.exists(paste0(results_path, "gg_sen_nonmon_data.RData"))) {
  gg_sen_nonmon_data <- data.table(
    dem_spend_adv = sen_X[sen_X[, "bottom_tail"] + sen_X[, "top_tail"] == 0, 1],
    d_dem_spend_adv_mean =
      apply(senate_boot_svm_cond_ame, 1 , mean, na.rm = TRUE),
    d_dem_spend_adv_q025 =
      apply(senate_boot_svm_cond_ame, 1 , quantile, probs = .025, na.rm = TRUE),
    d_dem_spend_adv_q975 =
      apply(senate_boot_svm_cond_ame, 1 , quantile, probs = .975, na.rm = TRUE)
  )
  save(gg_sen_nonmon_data,
    file = paste0(results_path, "gg_sen_nonmon_data.RData"))
} else {
  load(paste0(results_path, "gg_sen_nonmon_data.RData"))
}

if (!file.exists(paste0(results_path, "gg_senate_min_max_data_svm.RData"))) {
  available_seats <- senate_data[ , .(available_seats = .N), year]
  gg_senate_min_max_data_svm <- merge(
    senate_svm_counterfactuals[actual == 0 &
        dem_spend_adv == sen_counterfactual_vals[2],
      .(n_dems_max_dem = sum(boot_y > .5)),
      .(boot, year)],
    senate_svm_counterfactuals[actual == 0 &
        dem_spend_adv == sen_counterfactual_vals[1],
      .(n_dems_max_rep = sum(boot_y > .5)),
      .(boot, year)],
    by = c("boot", "year"))
  gg_senate_min_max_data_svm[, seat_swing := n_dems_max_dem - n_dems_max_rep]
  gg_senate_min_max_data_svm <- merge(gg_senate_min_max_data_svm,
    unincluded_senate_seat_totals[, .(year, D, R)], by = "year")
  gg_senate_min_max_data_svm <- merge(gg_senate_min_max_data_svm,
    available_seats, by = "year")
  gg_senate_min_max_data_svm[,
    n_reps_max_dem := available_seats - n_dems_max_dem]
  gg_senate_min_max_data_svm[,
    n_reps_max_rep := available_seats - n_dems_max_rep]
  save(gg_senate_min_max_data_svm,
    file = paste0(results_path, "gg_senate_min_max_data_svm.RData"))
} else {
  load(paste0(results_path, "gg_senate_min_max_data_svm.RData"))
}

if (!file.exists(paste0(results_path, "gg_senate_zero_money_data_svm.RData"))) {
  available_seats <- senate_data[ , .(available_seats = .N), year]
  gg_senate_zero_money_data_svm <- merge(
    senate_svm_counterfactuals[actual == 1,
      .(n_dems_actual = sum(boot_y > .5)),
      .(boot, year)],
    senate_svm_counterfactuals[actual == 0 & dem_spend_adv == 0,
      .(n_dems_zero = sum(boot_y > .5)),
      .(boot, year)],
    by = c("boot", "year"))
  gg_senate_zero_money_data_svm[, seat_swing := n_dems_actual - n_dems_zero]
  gg_senate_zero_money_data_svm <- merge(gg_senate_zero_money_data_svm,
    unincluded_senate_seat_totals[, .(year, D, R)], by = "year")
  gg_senate_zero_money_data_svm <- merge(gg_senate_zero_money_data_svm,
    available_seats, by = "year")
  gg_senate_zero_money_data_svm[,
    n_reps_actual := available_seats - n_dems_actual]
  gg_senate_zero_money_data_svm[,
    n_reps_zero := available_seats - n_dems_zero]
  save(gg_senate_zero_money_data_svm,
    file = paste0(results_path, "gg_senate_zero_money_data_svm.RData"))
} else {
  load(paste0(results_path, "gg_senate_zero_money_data_svm.RData"))
}
#     House SVM figure data ----
if (!file.exists(paste0(results_path, "gg_hou_nonmon_data.RData"))) {
  gg_hou_nonmon_data <- data.table(
    dem_spend_adv = hou_X[hou_X[, "bottom_tail"] + hou_X[, "top_tail"] == 0, 1],
    d_dem_spend_adv_mean =
      apply(house_boot_svm_cond_ame, 1 , mean, na.rm = TRUE),
    d_dem_spend_adv_q025 =
      apply(house_boot_svm_cond_ame, 1 , quantile, probs = .025, na.rm = TRUE),
    d_dem_spend_adv_q975 =
      apply(house_boot_svm_cond_ame, 1 , quantile, probs = .975, na.rm = TRUE)
  )
  save(gg_hou_nonmon_data,
    file = paste0(results_path, "gg_hou_nonmon_data.RData"))
} else {
  load(paste0(results_path, "gg_hou_nonmon_data.RData"))
}

if (!file.exists(paste0(results_path, "gg_house_min_max_data_svm.RData"))) {
  available_seats <- house_data[ , .(available_seats = .N), year]
  gg_house_min_max_data_svm <- merge(
    house_svm_counterfactuals[actual == 0 &
        dem_spend_adv == hou_counterfactual_vals[2],
      .(n_dems_max_dem = sum(boot_y > .5)),
      .(boot, year)],
    house_svm_counterfactuals[actual == 0 &
        dem_spend_adv == hou_counterfactual_vals[1],
      .(n_dems_max_rep = sum(boot_y > .5)),
      .(boot, year)],
    by = c("boot", "year"))
  gg_house_min_max_data_svm[, seat_swing := n_dems_max_dem - n_dems_max_rep]
  gg_house_min_max_data_svm <- merge(gg_house_min_max_data_svm,
    unincluded_house_seat_totals[, .(year, D, R)], by = "year")
  gg_house_min_max_data_svm <- merge(gg_house_min_max_data_svm,
    available_seats, by = "year")
  gg_house_min_max_data_svm[,
    n_reps_max_dem := available_seats - n_dems_max_dem]
  gg_house_min_max_data_svm[,
    n_reps_max_rep := available_seats - n_dems_max_rep]
  save(gg_house_min_max_data_svm,
    file = paste0(results_path, "gg_house_min_max_data_svm.RData"))
} else {
  load(paste0(results_path, "gg_house_min_max_data_svm.RData"))
}

if (!file.exists(paste0(results_path, "gg_house_zero_money_data_svm.RData"))) {
  available_seats <- house_data[ , .(available_seats = .N), year]
  gg_house_zero_money_data_svm <- merge(
    house_svm_counterfactuals[actual == 1,
      .(n_dems_actual = sum(boot_y > .5)),
      .(boot, year)],
    house_svm_counterfactuals[actual == 0 & dem_spend_adv == 0,
      .(n_dems_zero = sum(boot_y > .5)),
      .(boot, year)],
    by = c("boot", "year"))
  gg_house_zero_money_data_svm[, seat_swing := n_dems_actual - n_dems_zero]
  gg_house_zero_money_data_svm <- merge(gg_house_zero_money_data_svm,
    unincluded_house_seat_totals[, .(year, D, R)], by = "year")
  gg_house_zero_money_data_svm <- merge(gg_house_zero_money_data_svm,
    available_seats, by = "year")
  gg_house_zero_money_data_svm[,
    n_reps_actual := available_seats - n_dems_actual]
  gg_house_zero_money_data_svm[,
    n_reps_zero := available_seats - n_dems_zero]
  save(gg_house_zero_money_data_svm,
    file = paste0(results_path, "gg_house_zero_money_data_svm.RData"))
} else {
  load(paste0(results_path, "gg_house_zero_money_data_svm.RData"))
}

#     Make SVM version of Figure 1 ----

fig_house_nonmonotonic_svm <-
  ggplot(gg_hou_nonmon_data, aes(dem_spend_adv, d_dem_spend_adv_mean)) +
  geom_hline(yintercept = 0, lty = 3, color = "gray") +
  # geom_point(data = house_svm_data[incl == TRUE],
  #   aes(dem_spend_adv, d_dem_spend_adv),
  #   alpha = .1, shape = 21, fill = "black", color = "white",
  #   stroke = 0) +
  geom_ribbon(aes(ymin = d_dem_spend_adv_q025, ymax = d_dem_spend_adv_q975),
    alpha = .35) +
  geom_line(color = "blue") +
  scale_x_continuous(breaks = seq(-1.5, 1.5, .5)) +
  scale_y_continuous(breaks = c(0, .03, .06, .09),
    labels = c("0%", "3%", "6%", "9%")) +
  xlab(" ") + ylab("Marginal Effect of Spending") + ggtitle("House") +
  coord_cartesian(expand = FALSE, xlim = c(-1.75, 1.75), ylim = c(-.015, .09))

fig_senate_nonmonotonic_svm <-
  ggplot(gg_sen_nonmon_data, aes(dem_spend_adv, d_dem_spend_adv_mean)) +
  geom_hline(yintercept = 0, lty = 3, color = "gray") +
  # geom_point(data = senate_svm_data[incl == TRUE],
  #   aes(dem_spend_adv, d_dem_spend_adv),
  #   alpha = .5, shape = 21, fill = "black", color = "white",
  #   stroke = 0) +
  geom_ribbon(aes(ymin = d_dem_spend_adv_q025, ymax = d_dem_spend_adv_q975),
    alpha = .35) +
  geom_line(color = "blue") +
  scale_x_continuous(breaks = seq(-10, 10, 5)) +
  scale_y_continuous(breaks = c(0, .005, .01),
    labels = c("0%", "0.5%", "1%")) +
  xlab(" ") + ylab(" ") + ggtitle("Senate") +
  coord_cartesian(expand = FALSE, xlim = c(-11, 11), ylim = c(-.015/9, .01))

cairo_pdf(paste0("figures/fig-nonmonotonic-svm-", Sys.Date(), ".pdf"),
  width = 6.5, height = 4)
h <- .125
h2 <- .03
suppressWarnings(
  cowplot::ggdraw() +
    cowplot::draw_plot(fig_house_nonmonotonic_svm + theme(
      plot.title = element_text(size = 12),
      text = element_text(family = "CM Roman")),
      x = 0, y = h,  height = 1 - 2 * h, width = .5) +
    cowplot::draw_plot(fig_senate_nonmonotonic_svm + theme(
      plot.title = element_text(size = 12),
      text = element_text(family = "CM Roman")),
      x = .5, y = h, height = 1 - 2 * h, width = .5 - h2) +
    cowplot::draw_label(
      paste0("Nonlinear Effect Replication with SVM"),
      size = 18, x = .5, y = 1 - h / 2, fontface = "bold",
      fontfamily = "CM Roman") +
    cowplot::draw_label(paste0(
      "Democratic Candidate's Spending Advantage\n(in millions of ",
      "inflation-adjusted US$)"),
      size = 12, x = .5, y = .75 * h,
      fontfamily = "CM Roman")
)
dev.off()

#     Make SVM version of Figure 2 ----
fig_senate_min_max_svm <- ggplot(
  gg_senate_min_max_data_svm,
  aes(y = year, group = year, fill = ..x..)) +
  geom_vline(xintercept = 50, lty = 3, color = "gray") +
  geom_density_ridges(
    aes(x = D + n_dems_max_dem),
    bandwidth = 1,
    scale = 1, size = 0.1, rel_min_height = 0.01,
    color = "gray92",
    #alpha = .25,
    fill = "#0276FD") +
  geom_density_ridges(
    aes(x = D + n_dems_max_rep),
    bandwidth = 1,
    scale = 1, size = 0.1, rel_min_height = 0.01,
    color = "gray92",
    alpha = .5,
    fill = "#ef2f4f") +
  scale_x_continuous(limits = c(25, 75), expand = c(0, 0),
    breaks = seq(30, 70, 10)) +
  scale_y_reverse(breaks=seq(2016, 1980, -4), expand = c(0, 0.2)) +
  xlab("Democratic Senate Seats") +
  theme_minimal() +
  theme(panel.border = element_rect(color = "gray92", fill = NA),
    legend.position = "none",
    plot.title = element_text(hjust = .5)) +
  ggtitle("") +
  ylab("")

fig_house_min_max_svm <- ggplot(
  gg_house_min_max_data_svm,
  aes(y = year, group = year, fill = ..x..)) +
  geom_vline(xintercept = 218, lty = 3, color = "gray") +
  geom_density_ridges(
    aes(x = D + n_dems_max_dem),
    bandwidth = 4.35,
    scale = 1, size = 0.1, rel_min_height = 0.01,
    color = "gray92",
    #alpha = .25,
    fill = "#0276FD") +
  geom_density_ridges(
    aes(x = D + n_dems_max_rep),
    bandwidth = 4.35,
    scale = 1, size = 0.1, rel_min_height = 0.01,
    color = "gray92",
    alpha = .5,
    fill = "#ef2f4f") +
  scale_x_continuous(limits = c(0, 435), expand = c(0, 0),
    breaks = c(70, 145, 218, 290, 360)) +
  scale_y_reverse(breaks=seq(2016, 1980, -4), expand = c(0, 0.2)) +
  xlab("Democratic House Seats") +
  theme_minimal() +
  theme(panel.border = element_rect(color = "gray92", fill = NA),
    legend.position = "none",
    plot.title = element_text(hjust = .5)) +
  ggtitle("") +
  ylab("")

cairo_pdf(paste0("figures/fig-min-max-svm-", Sys.Date(), ".pdf"),
  width = 6.5, height = 4)
h <- .15
h2 <- .03
suppressWarnings(
  cowplot::ggdraw() +
    cowplot::draw_plot(fig_house_min_max_svm + ggtitle("House") +
        xlab("") + theme(
          plot.title = element_text(size = 12),
          text = element_text(family = "CM Roman")) +
        scale_y_reverse(breaks = seq(2016, 1980, -6),
          minor_breaks = seq(2016, 1980, -2), expand = c(0, .3)) +
        theme(panel.grid.minor = element_line(size = .1),
          panel.grid.major.y = element_line(size = .1)),
      x = 0, y = .35 * h,  height = 1 - 1 * h, width = .5) +
    cowplot::draw_plot(fig_senate_min_max_svm + ggtitle("Senate") +
        xlab("") + theme(
          plot.title = element_text(size = 12),
          text = element_text(family = "CM Roman")) +
        scale_y_reverse(breaks = seq(2016, 1980, -6),
          minor_breaks = seq(2016, 1980, -2), expand = c(0, .3)) +
        theme(panel.grid.minor = element_line(size = .1),
          panel.grid.major.y = element_line(size = .1)),
      x = .5, y = .35 * h, height = 1 - 1 * h, width = .5 - h2) +
    cowplot::draw_label(
      paste0("Congressional Control Depends on Spending (SVM)"),
      size = 18, x = .5, y = 1 - .4 * h, fontface = "bold",
      fontfamily = "CM Roman") +
    # cowplot::draw_label(
    #   paste0(""),
    #   size = 14, x = .5, y = 1 - .85 * h, fontface = "bold",
    #   fontfamily = "CM Roman") +
    cowplot::draw_label(paste0(
      "Number of Seats Held by Democrats Under Hypothetical Spending Profiles"),
      size = 12, x = .5, y = .4 * h,
      fontfamily = "CM Roman")
)
dev.off()

#     Make SVM version of Figure 3 ----
senate_sig_years <- gg_senate_zero_money_data_svm[,
  .(year, V1 = jitter(n_dems_zero - n_dems_actual, factor = .001))
  ][, .(mean(V1), quantile(V1, .025), quantile(V1, .975)), year][
    which(sign(V2) == sign(V3) | V2 == 0 | V3 == 0), year]
gg_senate_zero_money_data_svm[,
  sig := as.numeric(year %in% senate_sig_years)]

fig_senate_zero_money_svm <- ggplot(
  gg_senate_zero_money_data_svm,
  aes(y = year, group = year, fill = ..x..)) +
  geom_vline(xintercept = 0, lty = 3, color = "gray") +
  geom_density_ridges(
    aes(x = n_dems_zero - n_dems_actual, fill = factor(sig)),
    bandwidth = 1,
    scale = 1, size = 0.1, rel_min_height = 0.05,
    color = "gray92") +
  scale_y_reverse(breaks=seq(2016, 1980, -4), expand = c(0, 0.2)) +
  xlab("Change in Democratic Senate Seats") +
  theme_minimal() +
  theme(panel.border = element_rect(color = "gray92", fill = NA),
    legend.position = "none",
    plot.title = element_text(hjust = .5)) +
  ggtitle("") +
  scale_fill_manual(values = c("0" = "gray72", "1" = "gray42")) +
  ylab("")

house_sig_years <- gg_house_zero_money_data_svm[,
  .(year, V1 = jitter(n_dems_zero - n_dems_actual, factor = .001))
  ][, .(mean(V1), quantile(V1, .025), quantile(V1, .975)), year][
    which(sign(V2) == sign(V3) | V2 == 0 | V3 == 0), year]
gg_house_zero_money_data_svm[,
  sig := as.numeric(year %in% house_sig_years)]

fig_house_zero_money_svm <- ggplot(
  gg_house_zero_money_data_svm,
  aes(y = year, group = year, fill = ..x..)) +
  geom_vline(xintercept = 0, lty = 3, color = "gray") +
  geom_density_ridges(
    aes(x = n_dems_zero - n_dems_actual, fill = factor(sig)),
    bandwidth = 4.35,
    scale = 1, size = 0.1, rel_min_height = 0.05,
    color = "gray92") +
  # scale_x_continuous(limits = c(50, 385), expand = c(0, 0),
  #   breaks = c(70, 145, 218, 290, 360)) +
  scale_y_reverse(breaks=seq(2016, 1980, -4), expand = c(0, 0.2)) +
  xlab("Change in Democratic House Seats") +
  theme_minimal() +
  theme(panel.border = element_rect(color = "gray92", fill = NA),
    legend.position = "none",
    plot.title = element_text(hjust = .5)) +
  ggtitle("") +
  scale_fill_manual(values = c("0" = "gray72", "1" = "gray42")) +
  ylab("")

cairo_pdf(paste0("figures/fig-zero-money-svm-", Sys.Date(), ".pdf"),
  width = 6.5, height = 4)
h <- .15
h2 <- .03
suppressWarnings(
  cowplot::ggdraw() +
    cowplot::draw_plot(fig_house_zero_money_svm + ggtitle("House") +
        xlab("") + theme(
          plot.title = element_text(size = 12),
          text = element_text(family = "CM Roman")) +
        scale_y_reverse(expand = c(0, .3)),
      x = 0, y = .35 * h,  height = 1 - 1.5 * h, width = .5) +
    cowplot::draw_plot(fig_senate_zero_money_svm + ggtitle("Senate") +
        xlab("") + theme(
          plot.title = element_text(size = 12),
          text = element_text(family = "CM Roman")) +
        scale_y_reverse(expand = c(0, .3)),
      x = .5, y = .35 * h, height = 1 - 1.5 * h, width = .5 - h2) +
    cowplot::draw_label(
      paste0("The Effect of Removing Money (SVM)"),
      size = 18, x = .5, y = 1 - .8 * h, fontface = "bold",
      fontfamily = "CM Roman") +
    cowplot::draw_label(paste0(
      "Additional Seats Held by Democrats under Zero Spending"),
      size = 12, x = .5, y = .25 * h,
      fontfamily = "CM Roman")
)
dev.off()


# Appendix I ----
#     Smaller Counterfactuals Figure
sc_summary <- rbindlist(senate_krls_counterfactuals)
sen_ggdat <- sc_summary[!type_or_quantile %in% c("zero", "actual"), .(
  q025 = quantile(n_dems, .025),
  q975 = quantile(n_dems, .975),
  med = median(n_dems)
), .(year, q = as.numeric(type_or_quantile))]
sen_ggdat <- merge(sen_ggdat, unincluded_senate_seat_totals[, .(year, D)],
  by = "year")
sen_ggdat[, `:=`(q025 = q025 + D, med = med + D, q975 = q975 + D)]
sen_ggdat[, real_dem_expenditure_advantage_w_outside := senate_data[,
  quantile(real_dem_expenditure_advantage_w_outside, q)]]
sen_ggdat[, group := ifelse(q <= .5, 1, 2)]
obs_n_dems_elected <- data.table(year = seq(1980, 2018, 2),
  obs_n_dems_elected = c(47, 46, 47, 55, 55, 56, 57, 48, 45, 45, 50, 49,
    45, 51, 59, 53, 55, 46, 48, 47))
sen_ggdat <- merge(sen_ggdat, obs_n_dems_elected, by = "year")

hc_summary <- rbindlist(house_krls_counterfactuals)
hou_ggdat <- hc_summary[!type_or_quantile %in% c("zero", "actual"), .(
  q025 = quantile(n_dems, .025), q975 = quantile(n_dems, .975),
  med = median(n_dems)), .(year, q = as.numeric(type_or_quantile))]
hou_ggdat <- merge(hou_ggdat, unincluded_house_seat_totals[, .(year, D)],
  by = "year")
hou_ggdat[, `:=`(q025 = q025 + D, med = med + D, q975 = q975 + D)]
hou_ggdat[, real_dem_expenditure_advantage_w_outside := house_data[,
  quantile(real_dem_expenditure_advantage_w_outside, q)]]
hou_ggdat[, group := ifelse(q <= .5, 1, 2)]
obs_n_dems_elected <- data.table(year = seq(1980, 2018, 2),
  obs_n_dems_elected = rev(c(235, 194, 188, 201, 179, 257, 233, 203, 206, 213,
    212, 208, 205, 259, 268, 260, 258, 253, 269, 243)))
hou_ggdat <- merge(hou_ggdat, obs_n_dems_elected, by = "year")
hou_ggdat[, real_dem_expenditure_advantage_w_outside := house_data[,
  quantile(real_dem_expenditure_advantage_w_outside, q, na.rm = TRUE)]]

fig_senate_appendix_counterfactual <- ggplot(sen_ggdat,
  aes(real_dem_expenditure_advantage_w_outside))+
  facet_wrap(~ year) +
  geom_ribbon(aes(ymin = q025, ymax = q975, group = group), color = "gray",
    alpha = .5) +
  geom_hline(yintercept = 50, color = "red", linetype = "dotted") +
  geom_hline(aes(yintercept = obs_n_dems_elected), color = "blue",
    linetype = "dashed") +
  xlab("") +
  ylab("") +
  geom_line(aes(y = med, group = group)) +
  ggtitle("Senate") +
  geom_rug(data = senate_data[
    real_dem_expenditure_advantage_w_outside >= sen_tail_lims[1] &
      real_dem_expenditure_advantage_w_outside <= sen_tail_lims[2]],
    aes(x = real_dem_expenditure_advantage_w_outside)) +
  scale_x_continuous(breaks = c(-10, 0, 10))

fig_house_appendix_counterfactual <- ggplot(hou_ggdat,
  aes(real_dem_expenditure_advantage_w_outside))+
  facet_wrap(~year)+
  geom_ribbon(aes(ymin = q025, ymax = q975, group = group),
    color = "gray", alpha = .5) +
  geom_hline(yintercept = 218, color = "red", linetype = "dotted")+
  geom_hline(aes(yintercept = obs_n_dems_elected), color = "blue",
    linetype = "dashed")+
  geom_line(aes(y = med, group = group)) +
  xlab("") +
  ylab("Number of Seats Held by Democrats") +
  ggtitle("House") +
  geom_rug(data = house_data[
    real_dem_expenditure_advantage_w_outside >= hou_tail_lims[1] &
      real_dem_expenditure_advantage_w_outside <= hou_tail_lims[2]],
    aes(x = real_dem_expenditure_advantage_w_outside), alpha = .25) +
  scale_y_continuous(breaks = c(218-100, 218, 218+100))


cairo_pdf(paste0("figures/fig-small-counterfactuals-", Sys.Date(), ".pdf"),
  width = 6.5, height = 8)
h <- .125
h2 <- .03
suppressWarnings(
  cowplot::ggdraw() +
    cowplot::draw_plot(fig_house_appendix_counterfactual + theme(
      legend.position = "none",
      plot.title = element_text(size = 12),
      text = element_text(family = "CM Roman")),
      x = 0, y = h,  height = 1 - 2 * h, width = .5) +
    cowplot::draw_plot(fig_senate_appendix_counterfactual + theme(
      plot.title = element_text(size = 12),
      text = element_text(family = "CM Roman")),
      x = .5, y = h, height = 1 - 2 * h, width = .5 - h2) +
    cowplot::draw_label(
      paste0("Congressional Control with Different Spending Advantages"),
      size = 16, x = .5, y = 1 - h / 2, fontface = "bold",
      fontfamily = "CM Roman") +
    cowplot::draw_label(paste0(
      "Democratic Candidate's Spending Advantage\n(in millions of ",
      "inflation-adjusted US$)"),
      size = 12, x = .5, y = .75 * h+.025,
      fontfamily = "CM Roman")
)
dev.off()